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Accomplishments/New  Findings: 

•  We  developed  a  nonlinear  dynamical  systems  approach  to  auditory  processing  using  gradient 
frequency  networks  of  neural  oscillators.  Oscillators  in  the  network  are  driven  by  external 
forcing  and  receive  input  from  other  oscillators.  Both  types  of  interaction  may  involve  linear 
and/or  nonlinear  coupling,  which  can  evolve  over  time  via  a  generalized  form  of  Hebbian 
plasticity. 

•  We  conducted  three  types  of  theoretical  analyses  to  understand  signal  processing,  pattern 
formation  and  adaptation  in  networks  of  neural  oscillators.  First,  we  analyzed  and  categorized 
the  driven  behaviors  of  canonical  oscillators  under  periodic  forcing  in  four  parameter  regimes. 

•  Second,  we  analyzed  coupled  oscillators  for  each  parameter  regime  identified  in  the  previous 
study.  These  studies  of  systems  of  coupled  oscillators  enabled  us  to  better  understand  and 
utilize  the  complex  pattern  forming  dynamics  of  GrFNNs. 

•  Third,  we  completed  several  studies  of  Hebbian  plasticity  in  gradient  frequency  neural 
networks.  We  expanded  a  single-frequency  network  with  linear  learning  rule  into  a  multi¬ 
frequency  network  with  a  nonlinear  learning  rule.  We  discovered  a  set  of  rich  dynamics  that 
are  not  observed  in  either  single-frequency  systems  or  systems  without  plasticity.  We  analyzed 
the  dynamics  of  two  simple  oscillator  systems  with  plastic  connections:  an  oscillator  with 
plastic  coupling  to  single  external  input  and  two  oscillators  connected  by  plastic  coupling. 

•  We  developed  a  computational  modeling  framework  for  gradient  frequency  neural  networks  in 
Matlab®:  the  GrFNN  (pronounced  griffin)  Toolbox.  The  toolbox  was  developed  collaboration 
with  the  Music  Dynamics  Laboratory  at  UConn  and  has  been  released  to  the  research 
community:  https://github.com/MusicDynamicsLab/GrFNNToolbox. 

•  We  subcontracted  with  Array  Fire  to  speed  up  our  computational  simulations  using  GPU 
acceleration. 

•  We  developed  a  C++  version  of  the  gradient  frequency  neural  network  code  that  functions  as 
an  application  program  interface  (API).  It  can  be  used  to  develop  end-user  applications  that 
run  on  CPUs,  GPU,  mobile  platforms  and  embedded  devices. 

•  Next,  we  developed  three  models  of  auditory  processing.  First,  we  developed  and  analyzed  a 
model  of  cochlear  dynamics  based  on  coupling  between  linear  mechanical  resonance  of  the 
basilar  membrane  and  critical  nonlinear  oscillations  of  outer  hair  cells.  We  obtained  analytical 
forms  for  auditory  tuning  curves  of  both  unidirectionally  and  bidirectionally  coupled  systems. 
The  tuning  curves  of  the  model  fit  auditory  nerve  tuning  curve  data  from  the  macaque  monkey 
well. 

•  We  developed  a  canonical  model  of  mode-locked  neural  oscillation  in  the  human  auditory 
brainstem.  We  showed  that  the  model  could  reproduce  frequency  following  responses  (FFRs) 
to  musical  intervals  recorded  noninvasively  in  humans. 

•  We  developed  a  model  of  cortical  phase  locking  to  auditory  rhythms.  We  successfully  tested 
the  model’s  predictions  in  behavioral  and  neurophysiological  experiments. 
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Summary: 

Many  systems  in  nature,  especially  biological  systems,  are  nonlinear  dynamical  systems. 
Auditory  perception  provides  a  stunning  example  of  just  how  powerful  such  systems  can  be. 
Humans  recognize  complex  acoustic  patterns  under  challenging  listening  conditions,  such  as  a 
voice  in  a  crowded  room  or  on  a  city  street.  We  quickly  learn  patterns  that  are  significant  to  us, 
such  as  the  sound  of  a  new  dog’s  bark  or  the  rhythm  of  a  samba.  We  segregate  sounds  from  one 
another,  so  that  we  can  attend  to  one  sound  while  suppressing  others.  These  and  other 
remarkable  capabilities  of  human  perception  arise  from  nonlinear  dynamical  processes  in  the 
auditory  system. 

We  have  developed  a  theoretical  framework  for  auditory  neural  processing  and  auditory 
perception.  Our  approach  models  the  auditory  system  as  a  dynamical  system  consisting  of 
oscillatory  neural  networks,  and  auditory  perception  as  stable  dynamic  patterns  formed  in  the 
networks  in  response  to  acoustic  signals.  Our  models  capture  neural  dynamics  using  canonical 
dynamical  systems.  Canonical  systems  are  not  detailed  neurophysiological  models;  they  are 
generic  models  that  capture  the  neurocomputational  properties  of  a  family  of  neurophysiological 
models  using  bifurcation  theory.  Such  models  apply  across  temporal  and  spatial  scales.  For 
example,  they  are  appropriate  for  describing  critical  nonlinear  oscillations  of  outer  hair  cells  in 
the  cochlea,  mode-locking  of  chopper  cells  to  sound  in  the  cochlear  nucleus,  and  entrainment  of 
cortical  oscillations  to  auditory  rhythms. 

We  have  made  significant  theoretical  advances  in  understanding  how  gradient  frequency 
neural  networks  (GrFNNs)  respond  to  acoustic  stimulation  and  how  such  networks  can  learn  and 
adapt  through  Hebbian  plasticity.  We  analyzed  and  categorized  the  driven  behaviors  of  canonical 
oscillators  under  periodic  forcing.  We  conducted  series  of  studies  of  systems  of  coupled 
oscillators  to  better  understand  and  utilize  the  complex  network  dynamics  of  GrFNNs.  We 
completed  a  comprehensive  study  of  Hebbian  plasticity  in  gradient  frequency  neural  networks. 
We  discovered  a  set  of  rich  dynamics  that  are  not  observed  in  either  single-frequency  systems  or 
systems  without  plasticity. 

We  created  a  computational  modeling  framework  in  Matlab,  called  the  GrFNN  Toolbox,  and 
we  have  made  our  models  publicly  available.  We  have  developed  C++  versions  of  the  gradient 
frequency  neural  network  code  that  can  be  used  to  develop  end-user  applications  to  run  on  CPUs, 
GPUs,  mobile  platforms  and  embedded  devices.  We  have  shown  that  our  models  can  predict  and 
explain  various  aspects  of  auditory  processing  and  perception  that  are  difficult  to  account  for  by 
more  traditional  models  based  on  linear  signal  processing  techniques. 

We  began  development  of  three  important  classes  of  auditory  models.  First,  we  developed 
and  analyzed  a  canonical  model  of  cochlear  dynamics  based  on  coupling  between  linear 
mechanical  resonance  of  the  basilar  membrane  and  critical  nonlinear  oscillations  in  the  organ  of 
Corti.  To  validate  this  model  we  compared  it  with  auditory-nerve  tuning-curve  data  from  the 
macaque  monkey.  Second,  we  developed  a  canonical  model  of  mode-locked  neural  oscillation  in 
the  human  auditory  brainstem.  The  model  successfully  predicted  complex  nonlinear  population 
responses  to  musical  intervals.  Third,  we  developed  a  model  of  cortical  phase  locking  to  auditory 
rhythms.  The  model  predicted  the  results  of  behavioral  and  neurophysiological  experiments. 

Our  models  are  consistent  with  neurophysiological  evidence  on  the  role  of  neural  oscillation 
at  various  levels  of  the  auditory  system,  and  they  explain  phenomena  that  other  computational 
models  fail  to  explain.  The  predictions  hold  for  an  entire  family  of  dynamical  systems,  not  only 
specific  physiological  systems.  Thus,  our  model  provides  a  theoretical  and  computational 
framework  for  the  next  generation  of  auditory  processing  devices. 
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1.  Neural  Processing  at  Multiple  Timescales 


Many  systems  in  nature,  especially  biological  systems,  are  nonlinear  dynamical  systems. 
This  report  describes  an  approach  to  auditory  modeling  using  canonical  dynamical  systems. 
Canonical  systems  are  not  detailed  neurophysiological  models;  they  are  generic  models  that 
capture  the  computational  properties  of  a  family  of  neurophysiological  models  using  bifurcation 
theory.  Our  model  captures  the  responses  of  oscillatory  neural  systems  driven  with  time  varying 
external  signals  (Large,  Almonte,  &  Velasco,  2010).  It  is  applicable  across  temporal  and  spatial 
scales,  from  pitches  to  pitch  sequences  to  rhythmic  patterns.  As  such,  our  models  are  appropriate 
for  describing  various  phenomena  in  the  auditory  system,  including  critical  nonlinear  oscillations 
of  outer  hair  cells  (e.g.,  Fredrickson-Hemsing,  Ji,  Bruinsma,  &  Bozovic,  2012),  mode-locking  of 
choppers  in  the  cochlear  nucleus  (e.g.,  Laudanski,  Coombes,  Palmer,  &  Sumner,  2010),  and 
entrainment  of  oscillations  in  auditory  cortex  (e.g.,  S.  Nozaradan,  Peretz,  Missal,  &  Mouraux, 
2011). 

This  report  is  organized  as  follows.  The  remainder  of  Section  1  provides  a  brief  introduction 
to  gradient  frequency  neural  networks,  our  canonical  model  for  auditory  dynamics.  Sections  2  - 
4  then  describe  a  comprehensive  set  of  studies  that  we  carried  out  to  understand  the  signal 
processing,  pattern  formation  and  adaptive  properties  of  GrFNNs.  Section  2  analyses  signal 
processing  properties.  Section  3  analyzes  the  properties  of  coupled  systems  of  oscillators, 
focusing  on  mode-locking,  which  enables  complex  patterns  to  form  in  the  networks  when 
stimulated  with  sound.  Section  4  analyses  plasticity  in  canonical  oscillators,  finding  a  set  of  rich 
dynamics  that  are  not  observed  in  either  single-frequency  systems  or  systems  without  plasticity. 
Finally,  in  Section  5  we  describe  three  modeling  projects  that  show  how  well  our  models  can 
capture  auditory  processing  at  various  temporal  and  spatial  scales. 

1.1.  Gradient  Frequency  Neural  Networks  (GrFNNs) 

Gradient  frequency  neural  oscillator  networks  (GrFNNs)  are  canonical  neural  oscillators 
arrayed  along  a  tonotopic  frequency  gradient,  like  filter  banks  (Patterson  et  al.,  1992).  Unlike 
filter  banks,  however,  GrFNNs  have  nonlinear  properties  that  are  appropriate  to  model  networks 
of  spiking  neurons  (e.g.,  Hodgkin  &  Huxley,  1952)  or  mean-field  descriptions  of  neural 
populations  (e.g.,  Wilson  &  Cowan,  1972)  near  a  Hopf  or  a  Bautin  bifurcation  (Hoppensteadt  & 
Izhikevich,  2001).  Dynamical  properties  associated  with  these  bifurcations  have  been  identified 
in  the  auditory  system  (e.g.,  Fredrickson-Hemsing  et  al.,  2012;  Laudanski  et  al.,  2010). 

Our  canonical  model  (Eq.  1)  captures  the  dynamical  properties  of  a  network  of  neural 
oscillators  driven  with  an  acoustic  signal,  x(t),  and  synaptic  connections,  aj,  between  all  possible 


pairs  of  oscillators. 


T%Zi  —  I  a  +  i27T  +  pi  +  \zi\ 


yc _ £* 

^  yl-v 


y/eZj  1  -  yfizi  1  -  y/ex(t)  1  -  y/eZi  (1) 


where  zi  is  a  complex-valued  state  variable  that  represents  the  amplitude  and  phase  of  the  ith 
oscillator  in  the  network,  and  the  overdot  denotes  time  derivative,  and  the  roman  i  is  the 
imaginary  unit.  The  time  constant,  n,  determines  the  natural  frequency,  r,  =  1/fi.  The  right-hand 
side  of  Equation  1  consists  of  intrinsic  terms  (all  before  the  summation)  and  input  terms  from 
coupling  within  the  network  and  external  drive  xft)  (all  after  the  summation).  Depending  on  the 
parameter  values,  a,  fi,  fh  and  e,  a  canonical  oscillator  exhibits  one  of  several  distinct  intrinsic 
behaviors  available  near  a  Hopf  bifurcation  or  a  Bautin  (a.k.a.  double  limit  cycle)  bifurcation. 
Stability  analysis  shows  that  there  are  four  parameter  regimes  for  canonical  oscillators,  each  with 
a  distinct  set  of  autonomous  and  driven  behaviors  (Kim  &  Large,  2015)  (see  Section  2). 
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One  key  difference  between  W  (b) 

the  response  of  nonlinear 
oscillators  and  the  response  of 
linear  bandpass  filters  is  that 
individual  oscillations  mode -lock 
to  periodic  stimuli.  Mode-locking 
means  that  k  cycles  of  an 
oscillation  lock  to  m  cycles  of  the 
stimulus,  where  k  and  m  are 
integers,  as  shown  in  Fig  1.1b  ( k 
=  m  =  1  is  called  phase-locking). 

Frequency  (of  the  stimulus 
relative  to  the  oscillation)  and 
coupling  strength  determine  the 

locking  mode  (Fig  1.1b).  Mode-locking  responses  in  an  active  oscillatory  network  produce 
frequencies  that  are  not  present  in  its  stimulus:  harmonics,  subharmonics,  integer  ratios,  and 
combinations  of  stimulus  frequencies  (Fig  1.1b).  This  predicts  certain  types  of  network 
interactions  could  carry  out  complex  computations  needed  to  explain  the  perception  of  structure, 
for  example  the  perception  of  pitch  at  event  timescales  (Meddis  &  OMard,  2006)  and  the 
perception  of  pulse  and  meter  at  rhythmic  timescales  (Large,  Herrera,  &  Velasco,  2015)  (see 
Section  3). 

A  second  key  aspect  of  our  canonical  oscillator  networks  is  that  they  can  learn  and  adapt 
via  Hebbian  plasticity  (Hoppensteadt  &  Izhikevich,  1996).  Equation  2,  a  Hebbian  learning  rule, 
determines  the  dynamics  of  plastic  connections  that  alter  their  amplitude  and  phase  depending  on 
the  amplitude  and  frequency  relationship  between  its  source  and  target  (see  Section  4). 


Figure  1.1.  (a)  Intrinsic  oscillatory  dynamics  can  arise  through  a  Hopf  bifurcation, 
where  a  =  0  is  the  critical  point,  between  damped  (left)  and  spontaneous  (right) 
oscillation,  (b)  Within  each  resonance  region  in  the  frequency-coupling  parameter 
space  (c:  coupling  strength,  f:  oscillator  intrisic  frequency,  f0:  input  frequency)  the 
oscillator  mode-locks  to  the  input  at  the  k:m  ratio  shown.  Insets  show  example 
inputs  and  traces  produced  by  the  model. 
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In  Eq.  2  aj  is  a  complex  variable  representing  the  amplitude  and  phase  of  the  plastic  connection 
from  the  jth  oscillator  to  the  ith  oscillator.  The  parameters  X,  jui,  /m,  k,  and  ec  determine  the 
dynamics  of  the  plasticity.  A  simulation  of  the  model  defined  by  Equations  1  and  2  shows  that, 

unlike  linear  filters,  canonical  oscillators 
resonate  not  only  to  frequencies  present 
in  the  external  signal  but  also  to  their 
harmonics,  subharmonics,  and  nonlinear 
combinations  (Fig  1.2b).  Plastic 
connections  grow  strong  between 
oscillators  if  they  exhibit  certain  simple 
frequency  relationships  (Fig  1.2c).  In 
the  example,  a  “missing  fundamental” 
stimulus  is  input  to  the  network  (200  ad 
300  Hz),  and  the  network  responds  with 
the  missing  fundamental  frequency  (100 
Hz;  Fig  1.2b)  as  well  as  other  related 
frequencies.  The  connection  matrix 
learns  amplitude  and  phase  relationships 
of  the  components  (Fig  1.2c  shows  only  amplitude,  i.e.  I  aj  I)  and  this  information  can  be  used  to 
“bind”  the  simultaneous  components  into  an  integrated  percept. 


External  signal 

Source  oscillator  natural  frequency  (Hz) 

Figure  1.2.  A  GrFNN  model  with  external  forcing  and  plastic  internal 
connections.  (A)  Schematic  of  the  model  structure.  (B)  Oscillator 
amplitudes  and  (C)  plastic  connection  amplitudes  after  stimulation  with  an 
external  signal  containing  200  and  300  Hz  sinusoids.  Note  the  harmonic 
structure  of  the  response,  including  the  100  Hz  “missing  fundamental”. 
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Summary:  The  model  described  by  Equations  1  and  2  is  a  generic  population-level  model 
(Hoppensteadt  &  Izhikevich,  1997)  that  captures  the  fundamental  dynamics  observed  in  neuron- 
level  models  (Brunei,  2000;  Stefanescu  &  Jirsa,  2008),  but  is  amenable  to  theoretical  and 
computational  analysis  (Aronson,  Ermentrout,  &  Kopell,  1990).  The  model  is  invariant  over 
temporal  and  spatial  scale,  and  in  Section  5  we  use  it  to  model  active  cochlear  resonance  to 
sound  (Lerud,  Kim,  Almonte,  Carney,  &  Large,  under  revision),  brainstem  phase-locking  to 
pitch  combinations  (Lerud,  Almonte,  Kim,  &  Large,  2014),  and  cortical  entrainment  to  rhythmic 
patterns  (e.g.,  Large  et  ah,  2015).  First,  however,  we  describe  the  theoretical  analyses  that  have 
enabled  us  to  create  such  models. 

2.  Signal  Processing  by  Neural  Oscillators 

Despite  its  simple  mathematical  form,  the  canonical  model  for  gradient  frequency  neural 
networks  is  still  difficult  to  analyze  in  its  entirety  because  its  dynamics  is  determined  by  complex 
interactions  among  multiple  network  components.  Oscillators  in  the  network  are  driven  by 
external  forcing  and  at  the  same  time  receive  input  from  other  oscillators,  and  both  types  of 
interaction  may  involve  linear  and/or  nonlinear  coupling  which  can  evolve  over  time  via  a 
generalized  form  of  Hebbian  plasticity  (Large,  2011;  Hoppensteadt  and  Izhikevich,  1996).  Our 
approach  is  to  analyze  individual  components  of  the  network  separately  and  attempt  to 
understand  its  overall  dynamics  as  a  combination  of  its  component  dynamics.  To  begin,  we 
analyzed  and  categorized  the  driven  behaviors  of  canonical  oscillators  under  periodic  forcing. 

Consider  the  following  differential  equation  describing  an  oscillator  in  the  canonical 
model  (or  simply,  a  canonical  oscillator)  driven  by  sinusoidal  forcing  of  fixed  frequency,  coo,  and 
amplitude,  F: 

z  =  z(a  +  iu  +  Pi\z\2  +  +  FeiuJot 

V  i-eM2/  (3) 


where  co  =  2 nf  is  the  radian  natural  frequency.  To  understand  the  response  of  a  gradient 
frequency  network,  we  focus  on  how  the  driven  state  of  an  oscillator  changes  as  a  function  of  its 
natural  frequency. 

The  autonomous  behavior  of  the  oscillator  (i.e.,  when  F  =  0)  is  readily  seen  when  it  is 
brought  to  polar  coordinates  using  z  =  rel<p.  Then,  the  amplitude  and  phase  dynamics  are 
described  by 
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Figure  2.1.  Autonomous  behavior  of  a  canonical  oscillator  in  different 
parameter  regimes.  Amplitude  vector  field  is  shown  for  (A)  a  critical  Hopf 
regime,  (B)  a  supercritical  Hopf  regime,  (C)  a  supercritical  double  limit  cycle 
regime,  and  (D)  a  subcritical  double  limit  cycle  regime.  Filled  circles  indicate 
stable  fixed  points  (attractors)  and  empty  circles  unstable  fixed  points  (repellers). 
Arrows  indicate  the  direction  of  trajectories  in  the  vector  field. 


arrow  indicates  (Fig  2.1  A).  An  oscillator  with  this  type  of  amplitude 


(f>  =  UJ. 

Depending  on  the  values 
of  a,  /h,  and  // 2,  the 
autonomous  amplitude 
vector  field  (the  first 
equation  above)  can  have 
one  of  four  distinct 
topologies.  When  f 
decreases  monotonically 
as  r  increases,  the  origin 
is  the  only  fixed  point 
which  is  stable  as  the 
vector  field  decays  to  zero 
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Figure  2.2.  Driven  behavior  of  a  critical  Hopf  oscillator.  (A)  Steady-state 
amplitude  and  relative  phase  as  a  function  of  frequency  difference  (a  =  0, 
=  -100,  /?2  =  0,  F  =  0.2),  with  vertical  dashed  lines  indicating  the 
frequency  differences  used  for  panels  B-E,  (B)  trajectories  attracted  to  a 
stable  node  in  the  (r,y/)  plane  starting  from  a  set  of  different  initial 
conditions  (£l/2n  =  0. 1),  (C)  trajectories  attracted  to  a  stable  spiral  (L2/2jz  = 
0.5),  (D)  relative  phase  plotted  over  time  for  a  trajectory  in  panel  B  (phase 
locking),  and  (E)  relative  phase  plotted  over  time  for  a  trajectory  in  panel 
C  (phase  locking).  Filled  circles  in  panels  B  and  C  indicate  stable  fixed 
points. 


while  oscillating  at  its  natural 
frequency.  A  representative 
parameter  regime  for  this  type 
is  the  critical  point  of  a 
supercritical  Hopf  bifurcation 
(a  =  0,  /?i  <  0).  (A  subcritical 
Hopf  bifurcation  occurs  when 
a  =  0  and  /?i  >  0.)  When  r 
increases  from  the  origin  and 
then  decreases  after  a  local 
maximum,  there  is  a  stable 
nonzero  fixed  point  while  the 
origin  is  rendered  unstable  (Fig 
2.  IB).  An  oscillator  of  this  type 
shows  spontaneous  oscillation 
at  the  amplitude  of  the  stable 
fixed  point  (unless  the  initial 
condition  is  zero).  The 
supercritical  branch  of  a 
supercritical  Hopf  bifurcation 
(a  >  0,  /h  <  0)  is  an  example. 


When  there  are  three  fixed  points  with  two  local  extrema,  two  of  the  fixed  points  are  stable, 
indicating  bistability  between  equilibrium  at  zero  and  spontaneous  oscillation  at  a  nonzero 
amplitude  (Fig  2. 1C).  As  the  local  maximum  in  the  vector  field  moves  below  the  r  axis  by,  say, 
decreasing  /?i,  the  two  nonzero  fixed  points  collide  and  vanish  (Fig  2. ID).  This  transition  is 
called  a  double  limit  cycle  (hereafter,  DLC)  bifurcation  since  it  involves  two  limit  cycles  (closed 
orbits)  in  the  (r,<p)  plane,  one  stable  and  the  other  unstable.  Thus,  we  call  the  regime  shown  in 
Fig  2. 1C  (a  <  0,  fh  >  0,  fh  <  0,  local  max  >  0)  supercritical  DLC  and  the  one  shown  in  Fig  2. ID 
(a  <  0,  /?i  >  0,  /?2  <  0,  local  max  <  0)  subcritical  DLC.  The  subcritical  DLC  regime  has  only  one 
stable  fixed  point  at  zero  but  is  different  from  the  critical  Hopf  regime  (Fig  2.1  A)  in  that  it  has  a 
local  maximum  in  the  vector  field. 

To  examine  how  a  canonical  oscillator  responds  to  external  forcing,  we  bring  Equation  3 
to  polar  coordinates,  again  using  z  =  rel<p,  and  express  its  dynamics  in  terms  of  the  relative  phase 
(//  =  <p  ~  coot  so  that  a  stable  fixed  point  in  (r,  ip)  indicates  a  phase-locked  state: 

„5 


r  =  ar  +  /?ir3  +  jzrppi  +  F  cos  ip 

ip  =  i}~  F  sin  ^ 


(4) 


where  £2  =  co  -  coo  is  the  frequency  difference  between  the  oscillator  and  the  input.  We  evaluate 
the  stability  of  fixed  point(s)  for  a  range  of  forcing  parameters  Q.  and  F  wide  enough  to 
encompass  all  possible  qualitatively  different  driven  behaviors  of  the  four  regimes  of  intrinsic 
parameters  introduced  above. 
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2.1  Critical  Hopf  Oscillators 

A  stability  analysis  shows  that  a  critical  Hopf  oscillator  phase-locks  to  sinusoidal  forcing  of  any 
frequency  and  amplitude.  For  fixed  forcing  amplitude,  the  steady-state  amplitude  r*  is  maximum 
when  the  forcing  frequency 
is  the  same  as  the  natural 
frequency  (i.e.,  Q  =  0),  for 
which  the  steady-state 
relative  phase  yt*  is  zero 
indicating  in-phase 

synchronization  (Fig  2.2A). 

As  the  natural  frequency 
and  the  forcing  frequency 
become  more  different,  r* 
decreases  monotonically 

and  approaches  zero  while 

_|_7T 

yt*  approaches31 2.  While  the 
fixed  point  ( r*,y /*)  remains 
stable  for  all  values  of  Q,  it 
changes  its  stability  type 
from  a  stable  node  to  a 
stable  spiral  as  |Q|  increases 
from  0.  It  is  clearly  seen  in 
the  (r,y/)  space  that  the  two 
attractors  have  distinct  local 
trajectories  (Fig  2.2B  and  C). 

The  way  r  and  y/  approach 

their  steady-state  values  in  time  (monotonic  vs.  oscillating  approach)  reflects  the  difference 
between  a  node  and  a  spiral  (only  the  relative  phase  is  shown  in  Fig  2.2D  and  E). 


Time 


Time 


Figure  2.3.  Driven  behavior  of  a  supercritical  Hopf  oscillator  under  weak 
forcing.  (A)  Steady-state  amplitude  and  relative  phase  as  a  function  of 
frequency  difference  (a  =  1,  Pi  =  -100,  P2  =  0,  F  =  0.02),  with  vertical 
dashed  lines  indicating  the  frequency  differences  used  for  panels  B-E,  (B) 
trajectories  attracted  to  a  stable  node  in  the  (r,\|/)  plane  (Q/2jt  =  0.02),  (C) 
trajectories  drawn  to  a  limit  cycle  (Q/2jt  =  0.04),  (D)  relative  phase  plotted 
over  time  for  a  trajectory  in  panel  B  (phase  locking),  and  (E)  relative  phase 
plotted  over  time  for  a  trajectory  in  panel  C  (phase  slip).  In  panels  B  and  C, 
filled  and  empty  circles  indicate  stable  and  unstable  fixed  points 
respectively,  and  red  lines  show  limit-cycle  orbits. 


2.2  Supercritical  Hopf  Oscillators 

For  a  supercritical  Hopf  oscillator  under  weak  forcing,  there  exist  three  steady-state  solutions  for 
small  frequency  differences,  two  of  which  are  a  saddle -node  pair,  and  just  one  unstable  solution 
for  large  frequency  differences  (Fig  2.3A).  As  the  frequency  difference  increases  from  zero  for  a 
fixed  forcing  amplitude,  the  saddle  and  node  are  lost  via  a  saddle-node  invariant-circle  (SNIC) 
bifurcation  (also  called  a  saddle-node  infinite-period  or  SNIPER  bifurcation),  which  leaves  a 
stable  (attracting)  limit-cycle  orbit  with  an  unstable  fixed  point  inside  (Fig  2.3A-C).  For  stronger 
forcing,  only  one  fixed  point  exists  for  all  values  of  frequency  difference,  but  it  changes  from  a 
stable  node  to  a  stable  spiral  then  to  an  unstable  spiral  as  |Q|  grows  from  zero  (Fig  2.4A).  Now 
the  phase-locking  boundary  is  at  the  transition  from  a  stable  spiral  to  an  unstable  spiral  (i.e.,  a 
Hopf  bifurcation  in  the  (r,yj)  space). 

Outside  the  phase-locking  range  for  strong  forcing,  the  driven  behavior  of  a  supercritical 
Hopf  oscillator  can  be  divided  into  two  categories.  Just  outside  the  Hopf  boundary,  the  relative 
phase  changes  over  time  but  is  bounded  and  does  not  traverse  the  full  2k  range  (Fig  2.4D),  which 
is  called  a  libration  (as  opposed  to  a  rotation).  When  averaged  over  time,  this  “phase-trapped” 
oscillation  has  the  same  mean  frequency  as  the  input  frequency,  so  it  can  be  described  as 
frequency  locking  without  phase  locking  (Hoppensteadt  and  Izhikevich,  1997).  As  |Q|  increases 
further,  the  relative  phase  starts  making  full  rotations  (Fig  2.4E).  At  this  point,  the  average 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


instantaneous  frequency  of  the  oscillator  is  different  from  the  input  frequency  and  approaches  the 
natural  frequency  as  |Q|  approaches  infinity.  The  existence  of  phase-trapped  libration  (thus, 
frequency  locking)  outside  the  phase-locking  boundary  is  a  distinct  feature  of  the  Hopf  boundary. 

The  SNIC  phase¬ 
locking  boundary  and  the  Hopf 
boundary  exist  only  for  weak 
and  strong  forcing  levels 
respectively,  but  the  two  types 
of  phase-locking  boundary 
coexist  for  a  small  range  of 
intermediate  forcing  level. 

This  region  of  the  forcing 
parameter  space  contains  a 
complicated  but  well-studied 
set  of  bifurcations  such  as  the 
Bogdanov-Takens  bifurcation 
and  the  cusp  point.  Also,  it  is 
worth  noting  that  the  same  set 
of  bifurcations  are  found  for 
other  periodically  driven 
nonlinear  oscillators  or 
populations  of  oscillators  such 
as  the  forced  van  der  Pol 
oscillator  (Holmes  &  Rand, 

1978)  and  the  forced  Kuramoto  model  (Childs  &  Strogatz,  2008).  However,  the  canonical  model 
analyzed  here,  with  its  simple  mathematical  form,  allows  closer  analytical  examination  than  is 
possible  for  more  complex  models. 

2.3  Supercritical  Double  Limit  Cycle  Oscillators 

In  the  interest  of  space,  here  we  will  present  a  summary  of  the  driven  behaviors  of  supercritical 
DLC  oscillators  and  subcritical  DLC  oscillators  (see  Kim  &  Large,  2015,  for  a  full  analysis). 
These  are  also  summarized  in  Table  1  and  Fig  2.5  below. 

When  driven  by  a  sinusoid,  a  supercritical  DLC  oscillator  shows  three  distinct  sets  of 
behaviors  depending  on  the  strength  of  the  forcing,  and  many  of  these  behaviors  involve 
bistability.  Under  weak  forcing,  a  supercritical  DLC  oscillator  has  two  stable  fixed  points  for 
small  frequency  differences  and  only  one  for  large  frequency  differences.  The  stable  fixed  point 
with  a  small  amplitude  exists  for  all  values  of  Q,  but  the  one  with  a  high  amplitude  (a  stable  node) 
is  lost  via  a  SNIC  bifurcation  and  leaves  a  limit-cycle  rotation  (phase  slip)  where  it  collides  with 
a  saddle  point.  For  intermediate  forcing  amplitudes,  the  saddle-node  pair  at  high  amplitudes  still 
exists  for  small  frequency  differences,  but  the  stable  fixed  point  at  low  amplitudes  exists  only  for 
large  frequency  differences.  There  is  a  range  of  intermediate  frequency  differences  for  which  no 
stable  fixed  point  exists  and  all  trajectories  are  attracted  to  a  limit-cycle  rotation  that  the  saddle- 
node  pair  leaves.  As  the  frequency  difference  increases,  the  fixed  point  inside  the  limit  cycle 
changes  from  an  unstable  node  to  an  unstable  spiral  and  eventually  to  a  stable  spiral  (i.e.  a 
subcritical  Hopf  bifurcation).  When  the  forcing  amplitude  is  further  increased,  only  one  fixed 
point  exists  for  any  value  of  frequency  difference.  The  driven  state  (r,y/)  of  a  strongly  forced 
supercritical  DLC  oscillator  goes  through  transitions  from  a  stable  node  (phase  locking),  a  stable 
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Figure  2.4.  Driven  behavior  of  a  supercritical  Hopf  oscillator  under  strong 
forcing.  (A)  Steady- state  amplitude  and  relative  phase  as  a  function  of 
frequency  difference  (a  =  1,  Pi  =  -100,  p2  =  0,  F  =  0.2),  with  vertical 
dashed  lines  indicating  the  frequency  differences  used  for  panels  B-E,  (B) 
trajectories  attracted  to  a  phase-trapped  libration  in  the  (r,\j/)  plane  (Q/27T  = 
0.5),  (C)  trajectories  attracted  to  a  rotation  (Q/2n  =  0.7),  (D)  relative  phase 
plotted  over  time  for  a  trajectory  in  panel  B  (phase-trapped  frequency 
locking  without  phase  locking),  and  (E)  relative  phase  plotted  over  time  for 
a  trajectory  in  panel  C  (phase  slip).  In  panels  B  and  C,  empty  circles 
indicate  unstable  fixed  points,  and  red  lines  show  limit-cycle  orbits. 
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spiral  (phase  locking),  a  libration  around  an  unstable  spiral  (frequency  locking  without  phase 
locking),  a  rotation  around  an  unstable  spiral  (phase  slip),  and  lastly  bistability  between  phase 
locking  on  a  stable  spiral  and  phase  slip  on  a  stable  limit  cycle.  So,  a  strongly  forced  supercritical 
DLC  oscillator  has  two  phase-locking  boundaries,  a  supercritical  Hopf  bifurcation  and  a 
subcritical  Hopf  bifurcation. 
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2.4  Subcritical  Double  Limit  Cycle  Oscillators 

Like  a  critical  Hopf  oscillator,  a  subcritical  DLC  oscillator  is  attracted  to  an  equilibrium  at  zero 
when  it  is  not  driven.  But  the  presence  of  a  local  maximum  in  the  amplitude  vector  field  makes 
its  driven  dynamics  more  varied  and  interesting  than  that  of  a  critical  Hopf  oscillator.  Like  a 
supercritical  DLC  oscillator,  a  subcritical  DLC  oscillator  exhibits  three  different  sets  of  driven 

behaviors  depending  on  the 
forcing  amplitude. 

For  weak  forcing,  it 
behaves  like  a  critical  Hopf 
oscillator,  with  its  driven  state 
attracted  to  a  stable  node  when 
|Q|  is  small  and  to  a  stable  spiral 
when  |Q|  is  large.  For 
intermediate  forcing  amplitudes, 
a  pair  of  fixed  points  appears  at 
high  amplitudes  and  they  are  lost 
via  a  saddle-node  bifurcation  at  a 
certain  frequency  difference,  but 
this  saddle-node  bifurcation  does 
not  leave  a  limit  cycle  like  a 
SNIC  bifurcation.  When  driven 
strongly,  a  subcritical  DLC 
oscillator  has  the  same  set  of 
fixed  points  as  a  supercritical 
DLC  oscillator — a  stable  node,  a 
stable  spiral,  an  unstable  spiral, 
and  a  stable  spiral  as  |Q| 
increases  from  zero.  A 
supercritical  Hopf  bifurcation 
occurs  at  the  first  phase-locking 
boundary,  where  a  stable  spiral 
turns  unstable  and  a  stable  limit 
cycle  grows  around  it.  However, 
the  limit  cycle  does  not  grow  into  a  rotation  that  encompasses  the  origin,  which  is  the  case  for  a 
supercritical  DLC  oscillator.  Instead,  it  shrinks  back  and  turns  into  a  stable  spiral  via  another 
supercritical  Hopf  bifurcation.  In  the  absence  of  a  SNIC  or  subcritical  Hopf  bifurcation,  a 
strongly  driven  subcritical  DLC  oscillator  shows  no  bistability  and,  since  the  only  non-locked 
behavior  is  a  libration,  it  either  phase-locks  or  frequency-locks  to  the  input  for  all  values  of  Q. 
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Figure  2.5.  Stability  regions  for  a  canonical  oscillator  under  sinusoidal 
forcing.  The  stability  of  driven  state  (r*,\|/*)  is  shown  as  a  function  of 
forcing  amplitude  and  frequency  difference  for  (A)  a  critical  Hopf 
oscillator  (a  =  0,  pi  =  -100,  (32=  0),  (B)  a  supercritical  Hopf  oscillator  (a 
=  1,  (3i  =  -100,  (32=  0),  (C)  a  supercritical  double  limit  cycle  oscillator  (a 
=  -1,  (3i  =4,  p2=  -1,  s  =  1),  and  (D)  a  subcritical  double  limit  cycle 
oscillator  (a  =  -1,  (3i  =  2.5,  p2  =  -1,  e  =  1).  The  color  indicates  the 
stability  type  of  a  stable  fixed  point  if  there  is  one  (purple  if  there  are 
two).  If  there  is  no  stable  fixed  point,  the  color  indicates  the  stability  of 
an  unstable  fixed  point.  Dashed  horizontal  lines  indicate  the  forcing 
amplitudes  used  for  Figs.  2. 2-2.4. 


2.5  Classification  of  Parameter  Regimes  by  Driven  Behavior 

We  can  classify  all  possible  parameter  settings  for  canonical  oscillators  into  four  regimes 
with  distinct  driven  behaviors  (Table  1  and  Fig  2.5).  Oscillators  with  an  autonomous  amplitude 
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vector  field  that  monotonically  decreases  from  zero  with  no  local  extremum  (Fig  2.1  A)  have  the 
same  set  of  driven  behaviors  as  a  critical  Hopf  oscillator  (a  =  0,  /h  <  0,  (h  =  0).  Linear  oscillators 
(a  <  0,  /h  =  0,  [h=  0)  belong  to  this  category.  Oscillators  whose  amplitude  vector  fields  have  one 
local  maximum  (Fig  2. IB)  have  the  same  set  of  driven  behaviors  and  bifurcations  as  a 
supercritical  Hopf  oscillator  (a  >  0,  /?i  <  0,  [h  =  0).  Oscillators  with  a  <  0,  /h  >  0,  and  [h  <  0  are 
divided  into  three  groups  depending  on  whether  the  local  maximum  of  the  amplitude  vector  field 
is  above  zero  (Fig  2. 1C,  a  supercritical  DLC  oscillator),  is  below  zero  (Fig  2. ID,  a  subcritical 
DLC  oscillator),  or  does  not  exist  (with  the  same  driven  behaviors  as  a  critical  Hopf  oscillator). 

Table  1.  Classification  of  parameter  regimes  by  driven  behavior 


a 

a 

02 

Local  extremum" 

Discussed  as 

Bifurcations^ 

- 

0 

0 

None 

None 

0 

- 

0 

Critical  Hopf 

0 

0 

- 

- 

- 

0 

- 

0 

- 

0 

— 

— 

- 

+ 

- 

(No  max) 

+ 

- 

0 

One 

Supercritical  Hopf 

SNIC  (low  F); 

+ 

0 

- 

Super-Hopf  (high  F) 

0 

+ 

- 

+ 

- 

- 

+ 

+ 

- 

+ 

Two  (max  >  0) 

Supercritical  DLC 

SNIC  (low  F); 

SNIC,  Sub-Hopf  (mid  F); 

Super-Hopf,  Sub-Hopf  (high  F) 

+ 

Two  (max  <  0) 

Subcritical  DLC 

None  (low  F); 

SN  (mid  F); 

Super-Hopf,  Super-Hopf  (high  F) 

SNIC, 

saddle-node 

bifurcation  on  an 

invariant  circle;  Super-Hopf, 

supercritical  Hopf  bifurcation;  Sub-Hopf, 

subcritical  Hopf  bifurcation;  SN,  saddle-node  bifurcation. 

"The  number  of  local  extrema  in  the  autonomous  amplitude  vector  field. 
b  Bifurcations  at  phase-locking  boundaries. 


3.  Pattern  Formation,  Analysis  and  Recognition 

As  shown  in  Section  1,  we  use  dynamic  patterns  of  multi-frequency  oscillations  formed  in 
gradient  frequency  neural  networks  to  explain  and  model  auditory  neural  processing  and  auditory 
perception.  Here  we  present  the  results  of  a  series  of  analysis  for  systems  of  coupled  oscillators 
which  were  done  to  better  understand  and  utilize  the  complex  network  dynamics  of  GrFNNs. 

3.1  Stability  Analysis  of  Coupled  Oscillators:  Phase-Locking 

Here  we  analyze  the  following  system  of  two  canonical  oscillators  coupled  to  each  other: 

ii  =  Zi  (a  +  icui  +  /3i  |-2i  |2  + 

22  =  22  (a  4-  ic u2  -I-  di  M2  4- 

whose  polar  form  is 


Mdfili)  +  CZo 
1-ebi  P  )  +  CZ 2 


+  CZ\ 


(5) 
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h  =  an  +  fiirf  + 
r2  =  ar2  +  Pir%  + 

^  =  Vro  'sm^. 


+  cr2  cos  tj> 

+  CTi  COS  ^ 


(6) 


where  a  =  ne',p', yr  =  <pi  —  <p2  and  Q.  =  coi~  an.  We  present  the  analysis  of  two  coupled  oscillators 
for  each  parameter  regime  identified  in  Section  2. 

Critical  Hopf  oscillators  (a  =  0,  Pi  <  0,  Pi  =  0).  Eliminating  cos  t//*  from  the  steady-state 
amplitude  equations  reveals  that  there  are  only  fixed  points  with  symmetric  amplitudes  (i.e.  ri*  = 
ri).  Then,  we  can  reduce  Eq.  6  to  a  two-dimensional  system  with  the  phase  equation 
independent  of  amplitude  (r  =  n  =  n): 


f  =  r-3  +  cr  cos  tp 

ip  =  —  2csinV’- 


— tt/2 


0 


n/2 


-n/2 


0 

V1 


tv2 


-ti/2 


0 


tV2 


Figure  3.1.  Vector  field  of  the  relative  phase  of  two  coupled  oscillator  with 
symmetric  amplitudes  when  0  <  Q  <  2c  (left),  12  =  2c  (middle),  and  Q  >  2c  (right). 


The  steady-state  solution 
of  the  phase  equation 
above  indicates  that 
there  is  a  phase-locking 
boundary  at  |Q|  =  2c. 
There  is  a  node-saddle 
pair  when  |Q|  <  2c, 
which  vanishes  at  the 
boundary  (Fig.  3.1). 
When  |Q|  >  2c,  there  is 

no  i//*  (i.e.  no  phase-locking)  and  y/  rotates  in  either  positive  or  negative  direction.  Note  that  y/ 

TC  TC 

changes  slowly  (a  “bottleneck”)  near  ip  =  -  when  Q  >  2c  and  near  ip  =  —  -  when  Q  <  -2c. 

The  steady-state  solution  of  the  amplitude  equation  shows  that  nonzero  steady-state 

amplitude  exists  for  |Q|  <  2c  and  it  drops  to 
zero  at  the  boundary  (Fig.  3.2).  The 
Jacobian  matrix  and  its  trace  and 
determinant  suggest  that  the  fixed  point 
always  lies  on  the  border  between  stable 
nodes  and  stable  spirals  (it  is  a  star  node 
with  both  eigenvalues  at  -2c).  Zero  solution 
ri  =  r2=  0  is  also  a  possibility  for  coupled 
critical  oscillators,  regardless  of  |Q|  being 
greater  or  smaller  than  2c.  A  perturbation 
analysis  shows  that  the  zero  solution  is 
unstable  for  |Q|  <  2c,  and  stable  for  |Q|  >  2c. 

Supercritical  Hopf  oscillators  (a  > 
0,  pi  <  0,  pi  =  0).  These  are  two  possible 
solutions  for  two  coupled  supercritical 
oscillators:  solutions  with  symmetric 


-  Stable  node 
Stable  spiral 
Unstable  node 
Unstable  spiral 

-  Saddle  point 


Figure  3.2.  Steady-state  amplitude  and  relative  phase  of 
two  coupled  critical  Hopf  oscillators  (a  =  0,  (3i  =  -100,  c  = 
1).  Though  the  fixed  points  are  colored  orange,  they  are  all 
star  nodes. 
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0.15 


Figure  3.3.  Steady-state  amplitudes  (both  symmetric  and 
asymmetric)  and  relative  phase  of  two  coupled 
supercritical  Hopf  oscillators  (a  =  1,  pi  =  -100,  c  =  1). 
The  thick  lines  indicate  fixed  points  with  symmetric 
amplitudes  and  the  thin  lines  those  with  asymmetric 
amplitudes. 


amplitudes  and  those  with  asymmetric 
amplitudes.  The  symmetric  solution  is 


r  = 


\ 


a±\/c2-^ 

ft 


Examination  of  the  Jacobian  matrix 
indicates  the  first  solution  (the  bigger  r *)  is 
always  a  stable  node  and  the  second  one,  if 
it  exists,  is  a  saddle  point.  The  stable  node  is 
lost  via  a  saddle-node  bifurcation  at  the 
locking  boundary,  which  is  predictable  from 
the  observation  that  the  relative  phase 
rotates  right  outside  the  boundary  (Fig.  3.3, 
bold  lines).  The  asymmetric  solutions  are 


r,  = 


a 

'2ft 


1±  \  1- 


4  c2 


a2  +  ft2 


and  the  other  amplitude  is  obtained  by 

r*2  r*2  _  a 

rl  +r2  ~~Jl. 


However,  an  examination  of  the  Jacobian  matrix  reveals  that  asymmetric  solutions  are  all  saddle 
points.  So,  the  only  stable  solutions  are  coupled  supercritical  Hopf  oscillators  are  symmetric 
solutions  (Fig.  3.3). 

Supercritical  DLC  oscillators  (a  <  0,  ft  >  0,  ft  <  0,  local  max  >  0).  Fig.  3.4  shows  the 
stability  of  fixed  points  for  supercritical  DLC  oscillators  with  two  different  levels  of  coupling 
strength.  Symmetric  solutions  (thick  lines)  exist  only  when  |Q|  <  2c.  Within  this  range,  there  is 
only  one  stable  nonzero  fixed  point  with  symmetric  amplitudes,  which  is  always  a  node.  They 
also  have  stable  fixed  points  with  asymmetric  amplitudes,  due  to  their  bistability.  When  coupling 
is  weak  enough,  stable  asymmetric  solutions  exist  for  all  values  of  f2.  For  stronger  coupling,  they 
exist  for  large  values  of  |Q|  only.  Three  distinct  behaviors  are  possible  for  Q  >  2c.  Zero  solution 
is  stable  so  that  both  oscillators  can  decay  to  zero  if  initial  amplitudes  are  small.  When  both 
initial  amplitudes  are  big  enough,  oscillators  can  stay  not  locked  to  each  other,  with  relative 
phase  y/  rotating.  If  one  oscillator  starts  from  small  initial  amplitude  and  the  other  from  big 
amplitude,  they  could  phase-lock  with  asymmetric  steady-state  amplitudes. 

Subcritical  DLC  oscillators  (a  <  0,  ft  >  0,  ft  <  0,  local  max  <  0).  Subcritical  DLC 


Q/2n  Q/2n 


Figure  3.4.  Steady-state  amplitude  and  relative  phase  of  two  coupled  supercritical  DLC  oscillators  (a  = 
— 1,  pi  =  4,  p2  =  -1)  with  weak  coupling  (c  =  0.2,  left)  and  strong  coupling  (c  =  1.2,  right).  Thick  lines 
represent  symmetric  solutions  and  thin  lines  asymmetric  ones.  The  dashed  lines  indicate  the  phase¬ 
locking  boundary  for  symmetric  behaviors  |Q|  =  2c. 
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oscillators  have  nonzero  symmetric  r*’s  for  a  subset  of  |Q|  <  2c.  When  coupling  is  weak  and 
satisfies  the  following  condition,  there  is  no  nonzero  r*  for  all  values  of  Q.  and  zero  solution  is  the 
only  stable  fixed  point: 


,2  <_  ^2/?2  —  0i  +  2 yj 02(02  —  0{)  _  ^ 


This  is  when  the  local  maximum  in  amplitude  vector  field  is  below  zero  even  for  0  =  0.  With 
stronger  coupling,  nonzero  symmetric  r*’  s  exist  for  the  following  range  of  Q,  which  is  a  subset 
of  |Q|  <  2c: 

mi<2.c 

As  shown  in  Fig.  3.5,  fixed  points  with  asymmetric  amplitudes  exist  for  coupled  subcritical  DLC 
oscillators.  Most  of  them  are  unstable  saddle  points,  but  with  strong  coupling  there  is  a  narrow 
range  of  |Q|  just  outside  |Q|  =  2c  where  stable  spirals  exist. 


0/2n  Q/2n 


Figure  3.5.  Steady-state  amplitude  and  relative  phase  of  two  coupled  subcritical  DLC  oscillators  (a  =  —  1, 
(3i  =  2.5,  [L  =  -1)  with  weak  coupling  (c  =  0.3,  left)  and  strong  coupling  (c  =  1.2,  right). 


3.2  Stability  Analysis  of  Coupled  Oscillators:  Mode-Locking 

When  two  oscillators  have  natural  frequencies  whose  ratio  is  close  to  an  integer  ratio,  they  can 
mode-lock  to  each  other  via  resonant  monomials.  Two  canonical  oscillators  mode-locking  to 
each  other  can  be  described  by 

(  z\  =  zi  (a  +  icui  +  /?i|zi|2  +  | j  ^  cz^z"1  1 

\  z2  =  z2  (a  +  iu2  +  0i \z2 |2  +  Y 

for  which  we  assume  02  <  0  in  order  to  make  the  system  stable  given  arbitrarily  high  k  and  m. 
We  bring  the  system  to  polar  coordinates  using  z.i=  ne1<p‘  and  define  t//  =  m<p\  -  k<pi  and  Q  =  man  - 
ka> 2  to  get 

rx  -  qti  +  0i r2  +  crl^r™*1  cosxp 

<  r2  =  ar2  +  0i  +  £ k+™  cr™^-1  cosip 

ip  =  —  efc+,2  cr™_2r2_2(fcr2  +  mr^)  sin  ip. 


Assuming  symmetric  solutions,  the  system  can  be  a  two-dimensional  system  of  r  (=  n  =  n) 
and  y/,  for  which  the  steady-state  solutions  can  be  obtained  by  solving 
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a  +  for*2  +  Yz^z)  +(f^)  " 


2c2r*2(fc+m-2) 


Q/2n 


0/2tt 


To  obtain  asymmetric  solutions,  we 
must  go  back  to  the  original  three- 
dimensional  system.  We  do  not 
present  the  details  of  the  procedure 
since  equations  involved  in  getting 
asymmetric  solutions  are  too 
complicated  to  list  here.  Since  mode¬ 
locking  dynamics  do  not  change 
qualitatively  for  different  k’s  and  m’s, 
a  general  overview  of  both  symmetric 
and  asymmetric  behaviors  is  presented 
below  for  k  =  2,  m  =  1  only. 

As  shown  in  Fig.  3.6,  both 
supercritical  Hopf  oscillators  and  critical  Hopf  oscillators  have  only  symmetric  solutions.  For 
both,  stability  is  lost  through  a  saddle-node  bifurcation.  Supercritical  DLC  oscillators  have  both 
symmetric  and  asymmetric  stable  solutions  (Fig.  3.7,  left).  Note  that  for  asymmetric  solutions, 
which  are  stable  spirals, r2  is  bigger  thanri.  There  is  no  asymmetric  solution  with  ri  >  r2.  This  is 
because  the  system  is  asymmetric  with  respect  to  r\  and  n  due  to  unequal  k  =  2  and  m  =  1. 


Figure  3.6.  Steady-state  amplitude  and  relative  phase  of  two 
oscillators  coupled  via  resonant  monomials  for  2: 1  mode¬ 
locking:  supercritical  Hopf  oscillators  (a  =  1,  (3i  =  —  1,  (32  =  —  1,  s 
=  0.9,  c  =  1,  left)  and  critical  Hopf  oscillators  (a  =  0,  Pi  =  -1,  p2 
=  -1,  e  =  0.9,  c  =  1,  right). 


- Stable  node 

Stable  spiral 
Unstable  node 
Unstable  spiral 

- Saddle  point 

- - ^ - 

01  ‘  1  - — 1 - 1 - 1 - 1 - 1 —  1  1 

-1  -08  -06  -04  -0.2  0  0.2  04  0.6  0.8  1 


Q/2TT 


Figure  3.7.  Steady-state  amplitude  and  relative  phase  of  two  supercritical  DLC  oscillators  (a  =  —  1,  Pi  =  4, 
P2  =  — 1,  e  =  0.9,  c  =  1,  left)  and  two  subcritical  DLC  oscillators  (a  =  -1,  pi  =  2.5,  p2  =  — 1,  e  =  0.9,  c  =  1, 
right),  each  coupled  via  resonant  monomials  for  2: 1  mode-locking. 


Subcritical  DLC  oscillators  also  have  both  symmetric  and  asymmetric  solutions,  but  only  a  half 
of  symmetric  solutions  are  stable  (Fig.  3.7,  right). 

4.  Hebbian  Learning  and  Scene  Analysis 

Hoppensteadt  and  Izhikevich  (1996,  1997)  showed  that  a  weakly  connected  network  of  neural 
oscillators  of  identical  natural  frequencies  can  memorize  phase  differences  between  them  if 
synaptic  connections  change  in  time  according  to  complex  Hebbian  learning  rule: 

{zi  =  ( Pi  +  «*>)  Zi  -  Zi\zi\ 2  +  ^2 CijZj,  i  =  1, . . . , n 

j=1 

=  — ‘7 Cij  "I"  hjZiZj,  i  /  j 
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where  =  d/dz,  z  =  st  is  ‘slow’  time,  and  y  and  hj  are  positive  real  numbers.  Note  that  the  learning 
rule  has  only  a  linear  intrinsic  term  and  an  input  term.  When  z,i  and  zj  oscillate  at  the  same 
frequency,  the  input  term  kijZiZj  becomes  a  complex  constant  whose  phase  is  the  phase 
difference  between  zj  and  zj-  So,  the  phase  of  cy  eventually  comes  to  match  the  phase  difference 
between  the  oscillators  it  connects. 

We  expanded  Hoppensteadt  and  Izhikevich’s  single-frequency  network  with  linear 
learning  rule  into  a  multi-frequency  network  with  nonlinear  learning  rule.  We  show  that  this 
expansion  gives  rise  to  a  set  of  rich  dynamics  that  are  not  observed  in  either  single-frequency 
systems  or  systems  without  plasticity.  Here,  we  analyze  the  dynamics  of  two  simplest  oscillator 
systems  with  plastic  connections:  an  oscillator  with  plastic  coupling  to  single  external  input  and 
two  oscillators  connected  by  plastic  coupling. 

4.1  Plastic  Coupling  to  External  Forcing 

Consider  a  single  oscillator  z  that  is  driven  by  external  forcing  x  via  plastic  coupling  c. 
Introducing  fully  expanded  intrinsic  terms  to  both  the  oscillator  equation  and  the  learning  rule 
gives  us  the  following  system: 

i  =  z  (a  +  iw  +  di\z?  +  +  cx 

c  =  c  (A  +  |c|2  +  )  +KZX 

where  k  is  a  positive  real  number  representing  learning  rate. 

From  the  form  of  the  learning  rule, 
we  can  treat  the  connection  c  as  an  oscillator 
whose  natural  frequency  is  zero.  (Note  that 
there  is  no  imaginary  number  like  i co  in  the 
intrinsic  part  of  the  equation.)  This  means 
that  c  resonates  maximally  when  its  input 
kzx  has  constant  phase,  which  happens 
when  z  and  x  oscillate  at  the  same  frequency. 
When  z  and  x  are  not  oscillating  at  the  same 
Figure  4.1.  Waveform  of  sinusoidal  forcing  x,  a  critical  Hopf  frequency,  On  the  Other  hand,  C  is  driven  by 
oscillator  z,  and  plastic  coupling  c.  an  input  that  oscillates  at  the  difference 

frequency  and,  given  this  frequency  is  not  too  far  from  zero,  we  can  expect  that  c  would  phase- 
lock  to  this  oscillating  input  and  oscillate  at  the  difference  frequency  (Fig.  4.1). 

To  carry  out  the  analysis  further,  we  transform  the  system  to  polar  coordinates.  Assuming 
sinusoidal  forcing  of  constant  amplitude  and  frequency  (x  =  FelS,  d  =  coo)  and  defining  z  =  rel<p 
and  c  =  Aeld,  it  becomes 

r  —  ar  +  dir3  +  +  FA  cos (6  +  •&  —  </>) 

<\>  =  uj  +  EA  sin(0  +  •&  —  <!>) 

<  .  .  5 

A  =  A  A  +  Hi  A3  +  cos  (0  —  d  —  8) 

6  =  sin(4>  -  d  -  6). 

Since  the  only  angle  that  appears  on  the  right-hand  side  of  the  above  equations  is  8  +  d  -  cp,  we 
define  it  as  y/  and  get  the  following  three-dimensional  system: 

r  =  ar  +  dir3  +  +  FA  cos  xp 

<  A  =  A  A  +  n  i  A3  +  l'Ff2/Ai  T  nFr  cosxp 
=  -fl-  sin  xp 
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where  Q  =  co  -  coo.  Since  the  fully  expanded  system  is  difficult,  if  not  impossible,  to  solve,  let  us 
examine  the  following  truncated  version  (i.e.  we  assume  fh  =  /U2=  0): 

r  =  ar  +  fi\ r3  +  FA  cos  xp 
<  A  —  XA  +  +  nFrcosip 

ip  -  -Q-  F(Krr/^4  ^  sin0. 


Now  that  both  the  oscillator  and  the  learning  rule  have  multiple  parameter  regimes  with 
distinct  behaviors,  we  are  going  to  briefly  describe  each  combination  of  the  regimes  for  oscillator 
and  learning  rule.  First,  let  us  examine  the  system  consisting  of  a  critical  Hopf  oscillator  and 
linear  learning  rule — we  call  it  a  critical-linear  system.  Fig.  4.2  shows  r*.  A*,  if/*  and  steady-state 
9  of  a  critical-linear  system  as  functions  of  Q.  The  critical-linear  system  has  a  stable  fixed  point 

for  all  values  of  Q,  a  stable  node  for  small  IQI  and 
a  stable  spiral  for  large  IQI.  Notice  that  the 
oscillator’s  instantaneous  frequency  is  slower 
than  its  natural  frequency  when  the  oscillator’s 
natural  frequency  is  greater  than  the  input 

frequency,  and  the  other  way  around  when  the 
natural  frequency  is  smaller  than  the  input 

frequency.  See  also  that  the  connection’s 
instantaneous  frequency,  0 ,  is  always  between 
zero  and  Q  (marked  by  the  red  dashed  line).  Since 
9  =  0  —  0)0  when  the  system  is  in  a  steady  state, 
this  indicates  that  0  is  always  between  co  and  coo. 
We  also  observe  in  Fig.  4.2  that  as  IQI  increases, 
9  approaches  Q  and  hence  0  approaches  co.  So, 
when  IQI  is  large,  z  oscillates  near  its  natural  frequency  while  c  oscillates  with  small  amplitude  to 
cope  with  the  big  difference  between  z’s  instantaneous  frequency  and  the  input  frequency. 

For  a  critical-critical  system,  stable  fixed  points  exist  only  for  small  values  of  IQI  and 
they  appear  to  lie  on  the  border  between  stable  nodes  and  stable  spirals.  Outside  the  locking 
range,  both  the  oscillator  and  the  connection  slowly  decay  to  zero  with  their  amplitude 
fluctuating.  The  critical- supercritical  system  shown  in  Fig.  4.3  exhibits  interesting  dynamics. 
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Figure  4.2.  Steady-state  solutions  and  their  stability  for 
a  critical-linear  driven  system  (a  =  0,  fh  =  -1,  [h  =  0,  s  = 
0,  X  =  -1,  jlii  =  0,  in  -  0,  sc=  0,  k-  1,  F  -  1). 
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Figure  4.3.  Steady-state  solutions  and  their  stability  for 
a  critical-supercritical  driven  system  (a  =  0,  /h  =  -1,  Pi  = 
0,  £  =  0,  X  =  1,  JUl  =  — 1,  JU2  =  0,  £c  =  0,  k  =  1,  F  =  1). 


Figure  4.4.  Steady- state  solutions  and  their  stability  for 
a  supercritical- supercritical  driven  system  ( a  =  1,  = 

—  l,  [h  =  0,  £  =  0,  A  =  1,  JUl=  “I,  JU2=  0,  £c=  0,  k  =  1,  F  = 
1). 
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Contrary  to  the  critical-linear  system  shown  earlier,  0  converges  to  zero  (thus  <p  approaches  coo) 
as  IQI  increases  toward  infinity.  This  is  because  the  connection  has  nonzero  spontaneous 
amplitude  so  that  its  amplitude  cannot  be  lowered  indefinitely.  The  supercritical-linear  and 
supercritical-critical  systems  have  similar  overall  dynamics  to  the  critical-linear  system,  except 
that  the  oscillators’  amplitude  converges  to  their  nonzero  spontaneous  amplitude  as  IQI  increases. 
The  supercritical-supercritical  system  shown  in  Fig.  4.4  has  stable  fixed  points  for  only  small  IQI 
and  the  locking  boundary  is  a  saddle-node  bifurcation  point.  Outside  the  locking  range,  both  the 
oscillator  and  the  connection  fluctuate  near  their  spontaneous  amplitudes  and  intrinsic 
frequencies. 

4.2  Two  Oscillators  with  Plastic  Coupling 

The  dynamics  of  two  canonical  oscillators  connected  through  plastic  coupling  can  be  described 
by 

ii  =  zx  (a  +  kdi  +  ft |*i |2  +  )  +  c12z2 

Z 2  =  Z2  (<*  +  lu>2  +  Pi  \z2 12  +  +  c21^1 

Cl2  =  Cl2  (a  +  yUl  |ci2|2  +  fffcg  )  +  KZiZ2 
,  C21  =  C21  (A  +  Ml|c2l|2  +  +  K22^1- 

Defining  z.i  =  neUp‘  and  aj  =  Aijew«,  the  above  system  is  transformed  to 

h  =  arx  +  Pxr\  +  +  4i2r2  cos(0i2  +  fa  ~  <t>i) 

4>i=u\  +  sin(0i2  +  </>2  -  <t>i) 

r2  =  or 2  +  P\rl  +  +  42iri  cos(02i  +  <£i  -  <fo) 

<  4>2=u2  +  sin(6>2i  +<f> l  -  fa) 

Ai2  =  \Ai2  +  A*1^12  +  +  KT\r2  cos(^>i  -  <j>2  -  0i2) 

0\2  =  sin(0i  —  (\>2  —  012) 

A2i  =  XA2\  +  +  nv2v i  COS (02  —  01  —  #21) 

x  cc-ri21 

#21  =  !SA^'  Sin(^2  ~  4>1  ~  ^21  )• 

Note  that  there  are  only  two  distinct  arguments  for  trigonometric  functions  on  the  right-hand  side 
of  the  equations  above.  Defining  them  as  1//12  =  du-cpi+cpi  and  if/ 21  =  #21  -  cpi  +  cp\  turns  the  above 
eight-dimensional  system  into  a  six-dimensional  one: 

h  =  ari  +  Pir\  +  +  Al2r2  cos  ^12 

r2  =  ar2  +  ftrg  +  +  A2xri  cos  V21 

ii2  =  A4i2  4  Mi  A\2  +  +  «rir2  cos  ^12 

42i  =  A42i  +  Mi^li  +  +  Kr2ri  cos  ^21 

012  =  -Cl  -  sin V»i2  +  sin ^21 

ip2l=Cl-  ri(r^f8l)  sin  ^21  +  sin  ^12 

where  Q  =  cui  -  C02. 

From  the  definition  of  1//12  and  if/21,  we  can  see  that  012  =  cp1  —  cp2  and  d2i  =  02  —  0i 
when  the  system  is  in  a  steady  state  (i.e.  xp12  =  02i  =  0)-  In  other  words,  the  connections 
oscillate  at  the  instantaneous  frequency  difference  of  the  oscillators  they  connect.  This  also 
means  that  when  two  oscillators  have  different  instantaneous  frequencies,  two  connections 
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should  have  instantaneous  frequencies  of  the  same 
magnitude  but  in  opposite  directions  (i.e.  012  = 
—92i).  Another  property  that  can  be  expected 
from  the  definition  of  1//12  and  1//21  is  that  steady- 
state  On,  621  and  (p\  -q>2  should  be  neutrally  stable 
when  Q  =  0,  as  was  observed  for  driven  oscillators 
with  plastic  coupling  to  input. 

We  can  reduce  the  dimension  of  the 
system  further  by  assuming  that  two  oscillators 
and  two  connections  show  symmetric  behaviors. 
This  is  a  valid  as  well  as  useful  assumption,  since 
stable  fixed  points  are  always  symmetric  for  many 
parameter  regimes.  However,  it  does  not  capture 
any  of  asymmetric  behaviors  that  can  be  stable 
attractors  for  supercritical  DLC  oscillators.  Using  the  substitutions  r=r\=r2,A  =  An  =  A21  and  »// 
=  \jj\  2  =  -1//21,  we  get  the  following  three-dimensional  system: 
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Figure  4.5.  Symmetric  steady-state  solutions  and  their 
stability  for  a  critical-critical  coupled  system  ( a  -  0,  /?i  = 
-1,  /h  =  0,  £  =  0,  X  -  0,  JUl  =  -1,  JU2  =  0,  Sc=  0,  k  =  1). 


r  =  ar  +  /?ir3  +  +  Ar  cos  ^ 

A  =  XA  +  ViA3  +  f 
i>  = -SI-  sin^. 


2  COS  V’ 


As  we  did  above  for  oscillators  with  plastic  coupling  to  external  forcing,  we  will  briefly 
examine  different  possible  combinations  of  the  regimes  for  oscillator  and  learning  rule.  The 
critical-critical  system  (Fig.  4.5)  has  stable  nonzero  symmetric  solutions  for  only  small  values  of 
IQI.  The  symmetric  steady-state  behavior  of  a  critical-supercritical  system  does  not  look  much 
different  from  that  of  a  critical-critical  system.  It  too  has  nonzero  stable  fixed  points  for  small  IQI 
only.  But  note  that  since  the  learning  rule  is  in  a  supercritical  Hopf  regime,  A*  is  always  greater 
than  its  spontaneous  amplitude.  The  supercritical-linear  system  is  similar  to  the  driven  system  in 
the  same  regime — a  stable  node  for  small  IQI  and  a  stable  spiral  for  large  IQI.  A  difference 
between  the  two  systems  is  that  the  coupled  system  shown  here  is  unstable  and  blows  up  for 
small  values  of  IQI  when  k  >  [j\X.  Fully  expanded  intrinsic  damping  terms  with  negative  ^2  and  /U2 
solve  this  problem.  The  supercritical-critical  system  shown  in  Fig.  4.6  has  also  two  types  of 
stable  fixed  points  (nodes  and  spirals)  but  these  exist  for  small  values  of  IQI.  It  appears,  however, 

that  for  high  enough  k  all  unstable  fixed  points  in 
the  figure  turn  stable  and  become  stable  spirals. 
The  supercritical-supercritical  coupled  system, 
like  its  driven  counterpart,  has  stable  nodes  of 
amplitudes  greater  than  spontaneous  amplitudes 
for  small  values  of  IQI  only. 

Here  we  showed  the  analysis  of  symmetric 
solutions  only.  Numerical  simulations  suggest, 
however,  that  stable  steady- state  behaviors  are 
always  symmetric  for  these  systems  consisting  of 
truncated  canonical  models  of  oscillators  and 
connections.  So,  we  did  not  miss  any  stable 
solutions  above.  Simulations  also  show  that  stable 
asymmetric  solutions  exist  when  supercritical 


Figure  4.6.  Symmetric  steady-state  solutions  and  their 
stability  for  a  supercritical-critical  coupled  system  (a  = 
1,  fii  =  — 1,  fh  =  0,  s  =  0,  X  =  0,  jui  =  -1,  in  -  0,  sc  -  0,  k- 
0.8). 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


DLC  regimes  are  used,  but  analysis  of  the  kind  we  did  for  truncated  models  is  difficult  to  do  for 
asymmetric  behaviors  because  it  requires  solving  the  six-dimensional  system  with  nonzero  fh 
and  /J2. 

5.  Computational  Models  of  Physiological  and  Behavioral  Data 

As  described  in  the  preceding  sections,  we  have  developed  a  generic  population-level  model 
that  is  amenable  to  theoretical  analysis.  We  have  made  significant  theoretical  advances  in 
understanding  signal  processing,  pattern  formation  and  plasticity  in  networks  of  neural 
oscillators.  In  this  section  we  describe  now  we  have  applied  this  approach  to  model  active 
cochlear  resonance  to  sound  (Lerud  et  al.,  under  revision),  brainstem  mode-locking  to  pitch 
combinations  (Lerud  et  al.,  2014),  and  cortical  entrainment  to  rhythmic  patterns  (Large,  et  al., 
2015) . 

5.1  A  Canonical  Nonlinear  Cochlear  Model 

We  developed  and  analyzed  a  canonical  model  of  cochlear  dynamics  based  on  coupling 
between  linear  mechanical  resonance  of  the  basilar  membrane  and  critical  nonlinear  oscillations 
in  the  organ  of  Corti.  Through  dynamical  analysis,  we  obtain  analytical  forms  for  auditory  tuning 
curves  of  both  unidirectional  and  bidirectional  systems.  To  validate  this  model  we  compared  it 
with  auditory  nerve  (AN)  tuning  curve  data  from  the  macaque  monkey  (Joris  et  al.,  2011).  We 
found  that  the  model  tuning  curves  fit  the  macaque  data  with  good  accuracy.  A  fuller  description 
if  this  model  is  forthcoming  (Lerud,  Kim,  Almonte,  Carney  &  Large,  under  revision). 

Modern  models  of  the  cochlea  focus  on  the  nonlinear  oscillatory  responses  of  outer  hair 
cells  (Julicher,  2001;  Kern  &  Stoop,  2003).  Outer  hair  cells  poised  at  or  near  oscillatory  (Hopf) 
instability  (see  Section  2)  are  thought  to  be  responsible  for  the  cochlea's  extreme  sensitivity, 
excellent  frequency  selectivity,  and  amplitude  compression  (Camalet,  Duke,  Julicher,  &  Prost, 
1999;  Eguiluz,  Ospeck,  Choe,  Hudspeth,  &  Magnasco,  2000).  Models  of  OHC  nonlinearities 
consist  of  dynamical  equations  in  the  form  of  critical  oscillators  that  capture  generic  aspects  of 
nonlinear  resonance  (Fredrickson-Hemsing  et  al.,  2012;  Hoppensteadt  &  Izhikevich,  1997).  Such 
models  use  the  normal  (truncated)  form  of  the  Hopf  bifurcation  discussed  above  in  Section  2, 
(Equation  3).  Equation  4  shows  system  behavior  in  terms  of  amplitude  r  and  phase  angle  (f>.  A 
canonical  cochlear  model  in  which  each  segment  of  the  cochlea  is  represented  by  Equation  3 
(and  equivalently,  4)  can  account  for  a  nontrivial  subset  of  cochlear  dynamics  such  as  sharp 
mechanical  frequency  tuning,  exquisite  sensitivity,  and  a  large  dynamic  range  (Eguiluz  et  al., 
2000). 

Starting  with  normal  form  models  of  outer  hair  cell  nonlinearities  as  a  theoretical 
framework,  we  developed  an  extended  canonical  model  taking  into  account  linear  basilar 
membrane  dynamics,  critical  nonlinear  outer  hair  cell  dynamics,  and  the  coupling  between  the 
two.  In  the  first  model,  linear  basilar  membrane  oscillators  drive  critical  nonlinear  outer  hair  cell 
oscillators.  In  a  second  model,  bidirectional  coupling  was  introduced,  such  that  the  nonlinear 
elements  reciprocally  drive  the  linear  filtering  elements.  Both  models  can  be  solved  to  determine 
how  threshold  tuning  properties  depend  on  parameters,  and  both  models  produce  tuning  curves 
that  closely  match  responses  measured  in  the  macaque  AN  (Joris  et  al.,  2011).  In  addition,  our 
analysis  shows  that  the  bidirectionally  coupled  model  produces  intrinsic  oscillations,  such  that 
near  the  empirically  measured  threshold  there  exists  a  bifurcation  boundary  between 
nonsynchronized  and  synchronized  physiological  responses  (e.g., Johnson,  1980). 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


5.1.1  Unidirectional  Model 

We  used  pairs  of  coupled  oscillators  to  model  the  dynamics  of  cochlear  segments.  In  each 
pair,  one  oscillator  represents  BM  displacement  dynamics,  and  the  other  represents  organ  of 
Corti  (OC)  dynamics,  including  the  outer  hair  cells,  the  tectorial  membrane,  and  other  supporting 
structures.  Input  to  the  complex  drives  the  BM  oscillator,  which  is  intended  to  account  for  the 
dynamical  effects  of  the  cochlear  fluid  traveling  waves  that  drive  the  BM.  The  OC  energy  source 
stems  from  critical  oscillations  that  cause  the  organ  of  Corti  to  vibrate.  Thus,  the  model  exhibits 
both  BM  filtering  and  critical  oscillations  that  capture  the  amplification,  compression,  and 
frequency  selectivity  of  cochlear  processing  (Eguiluz  et  al.,  2000). 

The  natural  frequency  of  each  BM-OC  complex  is  set  to  correspond  to  the  best  frequency 
of  the  cochlear  segment  that  it  represents.  We  equate  the  state  of  the  OC  oscillator  with  the  signal 
that  is  transmitted  to  the  AN.  These  broad  considerations  lead  to  a  coupled  set  of  canonical 
oscillator  equations  for  modeling  a  BM-OC  complex: 

Zbm  —  Zbm  (<^b m  +  ^7T /)  +  Fel2n^ot 
Zoc  —  Zoc  (^OC  -b  12tt/  +  (P  +  iS)  |  Zqc\  )  T  C2±Zf}m 


The  state  variable  zbm  represents  the  dynamics  of  the  BM,  while  zoc  represents  the 
dynamics  of  the  OC,  including  the  nonlinearities  of  the  OHCs.  For  simplicity,  we  assume  a  linear 
BM.  This  leads  to  bandpass  filtering  behavior,  making  the  model  conceptually  similar  to  that  of 
(Julicher,  Andor,  &  Duke,  2001).  The  linear  damping  parameter,  ccbm  <  0,  is  determined  by 
fitting  tuning-curve  data.  For  the  OC  we  assume  critical  nonlinear  oscillation,  i.e.,  aox  =  0, 
resulting  in  optimal  amplification.  The  nonlinear  damping  parameter  (3  <  0  provides  amplitude 
compression  in  the  OC  and  is  also  determined  by  fitting  tuning-curve  data.  Finally,  the  parameter 
C21  governs  the  relative  strength  of  forcing  of  the  OC  by  the  BM  and  is  determined  by  fitting  the 
data  as  well. 

Because  the  model  is  described  in  terms  of  the  complex  state  variables  zbm  and  zoc,  it  can 
be  rewritten  in  polar  form,  giving  rise  to  amplitude  and  phase  equations: 


I'brri  —  OlbmTbrn  T  F  COs(2,7T fot  (ftbm') 

F 

(j)brn  =  27 rf  H - sin(27r/0f  -  (j)brn ) 

T'bm 

foe  C^ocT'oc  H”  fiv'oc  "I”  ^21  T'bm  COS ((frbm  ^oc) 

4>oc  =  2t rf  +  C2irbm  sm((j)brn  -  <f>oc). 
foe 

Given  the  threshold  amplitude  r*oc,  which  is  a  small  number  that  we  hold  constant  across  tuning 
curves,  the  formula  for  F  only  in  terms  of  model  parameters  is 


F 


+  2  ro‘caocP  +  r**(32  +  O2  +  2  Cl5r 


*2 

OC 


+  S2r 


1 

2 


F  is  normally  in  pascals  which  can  then  be  converted  to  the  stimulus  level  L  in  dB  SPL  by  L=20 


log  (F/po)  -  G,  where  po  =  20pPa  represents  the  reference  pressure,  and  G  the  gain  of  the  middle 
ear  filter  in  dB. 
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To  determine  tuning  curves  that  can  be  compared  with  auditory-nerve  data,  we  first  pass 
the  acoustic  stimulus  through  a  linear  filter  to  approximate  the  amplitude  and  phase  response  of 
the  middle  ear  (Bruce,  Sachs,  &  Young,  2003;  Zilany  &  Bruce,  2006).  The  middle-ear  filter  is  a 
simplified  form  of  that  of  Bruce  (2003).  Zilany  and  Bruce  developed  a  fifth-order  continuous¬ 
time  transfer  function  and  represented  it  as  a  fifth-order  digital  filter  using  a  bilinear 
transformation  for  a  sampling  frequency  of  500  kHz,  with  the  frequency  axis  pre-warped  to  give 
a  matching  frequency  response  at  1  kHz.  To  ensure  stability  of  the  digital  filter,  it  was 
implemented  in  a  second-order  system  form  by  cascading  digital  filters.  Once  the  MEF  and 
cochlear  BM-OC  oscillatory  complexes  are  defined,  the  resulting  waveform  is  provided  as  input. 

The  three  parameters,  auu  <  0,  <  0,  and  C21  >  0,  are  determined  using  a  search 

procedure  that  adjusts  model  parameters  to  obtain  a  sufficiently  close  match  to  the  data.  We  held 
the  threshold  r*oc  =  0.1  constant  and  fit  each  curve  individually.  The  results  of  the  parameter 
searches  are  shown  in  green  in  Figure  5.1  A,  B,  and  C,  for  low,  mid,  and  high  frequency  tuning 
curves  from  the  Joris  et  al.  data  set,  respectively.  The  average  root-mean-square  error  for  the  fits 
in  was  6.9219  dB.  It  was  noted  that  both  ccbm  and  /?  varied  within  a  single  order  of  magnitude 
over  all  tuning  curves,  and  cn  was  a  reliable  linear  function  of  CF. 

5.1.2  Bidirectional  Model 

A  more  realistic  configuration  of  the  BM-OC  oscillatory  complexes  is  to  use  bidirectional 
coupling  between  the  two  oscillators  rather  than  unidirectional  coupling  used  in  the  previous 
model.  Thus,  our  second  model  considers  the  effect  of  OC  dynamics  on  the  BM.  With 
bidirectional  coupling,  the  dynamics  of  an  oscillatory  complex  are  governed  by 

Zbm  =  zbm  ( abm  +  i2n f)  +  Fel27Tfot  C]_2^oc 

%oc  Zoc  (c^oc  +  127 r/  +  (/3  +  i5)|zoc| 2)  +  C21  zbm, 

with  C12  being  the  coupling  coefficient  of  the  OC  oscillator  feeding  back  to  the  BM  oscillator. 
Similarly  to  the  unidirectional  model,  we  can  get  a  closed-form  formula  for  forcing  amplitude  F 
expressed  as  a  function  of  threshold  amplitude  r*oc,  frequency  difference  Q,  and  model 
parameters: 

F  =  \j ~{^bnFlrn  +  COS  +  (ftr*m  +  C12^c  sin^*c)2 


where 


bm 


=  —\j  (ooc  +  /3r*2)2  +  (ft  +  <ir*2)2 , 

C21 


sin  ip: 

cos 


* 

OC 


r*c  (ft  +  5r%) 


C21^ 


and 


\J  1  —  sin2 


Stability  analysis  of  this  system  reveals  that  bidirectional  coupling  introduces  an  important 
change  in  the  dynamics  of  BM-OC  complexes.  With  unidirectional  coupling  and  with  a\m  set  to 
zero  or  below,  both  the  BM  and  OC  oscillators  decay  to  zero  when  not  driven  by  external  forcing. 
With  bidirectional  coupling,  however,  the  two  oscillators  provide  input  to  each  other  and  as  a 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


result  they  have  nonzero  steady-state  amplitudes  even  in  the  absence  of  external  forcing.  A 
consequence  of  having  nonzero  spontaneous  amplitude  is  that  a  BM-OC  complex  with 
bidirectional  coupling  may  not  phase-lock  to  external  forcing  if  its  natural  frequency  is  too 
different  from  forcing  frequency  or  if  forcing  is  not  strong  enough.  Figures  5.1D,E&F  show 
resonance  regions  or  “Arnold  tongues”  for  a  bidirectional  model  within  which  the  model  phase- 
locks  to  external  forcing.  Steady-state  solutions  (or  fixed  points)  may  exist  outside  the  resonance 
region,  but  they  are  not  stable  (i.e.,  the  model  is  not  attracted  to  them).  Typically,  steady-state 
amplitudes  r*oc  and  r* bm  are  unstable  when  they  are  smaller  than  the  spontaneous  amplitudes.  In 
comparison,  a  BM-OC  complex  with  unidirectional  coupling  and  ana  <  0  phase-locks  to 
external  forcing  of  any  frequency  and  amplitude. 

Due  to  the 
possibility  of 

unstable  solutions 
for  bidirectional 
models,  the  forcing 
amplitude  F  should 
be  examined  for  its 
stability  when  the 
threshold  amplitude 
r*oc  set  below  the 
spontaneous 
amplitude  of  zoc.  To 
compare  tuning  in 
this  model  to  the 
tuning  in  the 
unidirectional 
model,  we  choose 
the  coupling,  ci2, 
such  that 

spontaneous 
amplitude  is  just 
below  threshold 
amplitude,  r*oc  = 

0.1.  Figure  5.1 
D,E&F  show  the 

tuning  curves  (red)  lie  just  above  the  phase-locking  boundary,  and  are  similar  to  the  tuning 
curves  for  the  unidirectional  model.  The  bidirectional  model  did  not  provide  significantly  better 
fits  to  the  tuning-curve  data.  However,  the  bidirectional  model  makes  the  important  prediction 
that  cochlear  and  AN  phase-locking  will  be  observed  before  amplitude  of  firing  frequency 
increases  significantly.  This  prediction  matches  empirical  observations  (e.g.,  Johnson,  1980). 
Another  important  feature  of  this  bidirectional  model  is  spontaneous  oscillation.  It  is  well  known 
that  many  mammalian  cochleae  exhibit  spontaneous  otoacoustic  emissions,  and  this  aspect  of 
auditory  nonlinearity  cannot  be  predicted  with  a  unidirectionally  coupled  model. 

Summary:  Nonlinear  responses  to  acoustic  signals  arise  through  active  processes  in  the 
cochlea,  whose  exquisite  sensitivity  and  wide  dynamic  range  can  be  explained  by  critical 
nonlinear  oscillations  of  hair  cells.  We  studied  how  the  interaction  of  critical  nonlinearities  with 
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Figure  5.1.  Top:  Fits  of  the  unidirectional  model  to  low  (A),  mid  (B),  and  high  (C) 
frequency  AN  fibers  from  the  Joris  (2001)  data  set.  Data  is  in  blue,  model  fits  are  in 
green.  Bottom:  Resonance  regions  of  the  bidirectional  model  for  low  (D),  mid  (E),  and 
high  (F)  frequency  AN  fibers,  using  the  same  parameters.  Coupling  from  OC  to  BM, 
was  chosen  so  that  spontaneous  amplitude  was  slightly  below  threshold  amplitude, 
roc  =  0.1.  The  red  contours  show  threshold  amplitude.  The  BM-OC  model  phase- 
locks  to  external  forcing  in  the  parameter  regions  where  the  fixed  point  is  either  a 
stable  node  (red)  or  a  stable  spiral  (yellow).  Non-phase  locked  regions  (saddle  points) 
are  shown  in  blue. 
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the  basilar  membrane  and  other  organ  of  Corti  components  could  determine  tuning  properties  of 
the  mammalian  cochlea.  We  developed  a  canonical  model  in  which  the  dynamics  of  the  basilar 
membrane-organ  of  Corti  interaction  is  captured  using  pairs  of  coupled  oscillators  tuned  to  a 
gradient  of  natural  frequencies.  We  first  developed  a  minimal  model  in  which  a  linear  oscillator, 
representing  basilar  membrane  dynamics,  is  coupled  to  a  nonlinear  oscillator  poised  at  a  Hopf 
instability,  which  captures  the  nonlinear  responses  of  outer  hair  cells  and  related  organ  of  Corti 
components.  Parameters  were  determined  by  fitting  the  auditory-nerve  tuning  curves  of  macaque 
monkeys.  We  then  developed  a  more  sophisticated  model,  taking  into  account  bidirectional 
coupling.  We  found  that  the  unidirectionally  and  bidirectionally  coupled  models  account  equally 
well  for  threshold  tuning,  but  the  bidirectionally  coupled  model  also  exhibited  low  amplitude, 
spontaneous  oscillation,  providing  a  model  that  phase-locks  to  sound. 

5.2  Brainstem  Processing  of  Pitch  Combinations 

While  some  nonlinear  responses  arise  through  active  processes  in  the  cochlea,  others  arise 
in  neural  populations  of  the  cochlear  nucleus,  inferior  colliculus  and  higher  auditory  areas. 
Recently  mode-locking,  a  generalization  of  phase  locking  that  implies  an  intrinsically  nonlinear 
processing  of  sound,  has  been  observed  in  mammalian  auditory  brainstem  nuclei.  We  developed 
a  canonical  model  of  mode-locked  neural  oscillation  in  brainstem  that  predicts  the  complex 
nonlinear  population  responses  to  musical  intervals  that  have  been  observed  in  the  human 
brainstem  (for  complete  details,  see  Lerud  et  al.,  2014). 

In  central  auditory  circuits,  action  potentials  phase-lock  to  both  the  fine  time  structure 
and  the  temporal  envelope  modulations  of  auditory  stimuli  at  many  different  levels,  including 
cochlear  nucleus,  superior  olive,  inferior  colliculus,  thalamus,  and  Al  (Langner,  1992). 
Traditionally,  phase-locked  spiking  in  the  central  auditory  system  is  thought  to  represent  an 
essentially  passive  transmission  of  synchronized  basilar  membrane  motion.  An  alternative 
possibility  is  that  active  circuits  in  the  central  auditory  system  carry  synchronized  neural  activity 
forward.  If  this  is  the  case,  nonlinearities  observed  at  the  level  of  the  brainstem  might  also  arise 
due  to  mode-locking  (see  Section  3),  a  phenomenon  that  has  been  observed  in  the  auditory 
brainstem  (Langner,  1992),  and  is  physiologically  distinct  from  the  mechanical  compression  and 
half-wave  rectification  that  occurs  in  the  organ  of  Corti. 

Mode-locking  to  acoustic  signals  has  been  observed  in  guinea  pig  cochlear  nucleus 
chopper  and  onset  neurons  (Laudanski  et  al.,  2010),  and  mode-  locking  to  the  difference  tone  of 
two  dichotically  presented  stimulus  frequencies  has  been  observed  in  vivo  and  isolated  to  the 
inferior  colliculus  of  the  chinchilla  (Arnold  and  Burkard,  1998,  2000).  Mode-locked  spiking 
patterns  are  often  observed  in  vitro  under  DC  injection  (Brumberg  and  Gutkin,  2007),  and  active 
oscillations  have  been  observed  in  vivo  in  the  inferior  colliculus  of  the  chicken  (Schwarz  et  al., 
1993).  Such  observations  lead  to  the  possibility  that  the  nonlinear  responses  observed  in  the 
human  auditory  brainstem  might  arise,  in  part,  due  to  mode-locking  neurodynamics. 

We  modeled  nonlinear  responses  to  musical  intervals  that  have  been  measured  in  the 
human  auditory  brainstem  response  (Lee  et  al.,  2009,  see  Fig.  2).  In  that  study,  the  brainstem 
representation  of  the  musical  intervals  comprised  not  only  stimulus  frequencies,  but  also 
numerous  resonances  at  frequencies  that  were  not  physically  present  in  the  stimulus.  How  did 
these  frequencies  arise?  The  stimuli  were  the  intervals  major  sixth  (G  and  E,  “consonant”)  and 
minor  seventh  (F#  and  E,  “dissonant”)  which  have  fundamental  frequency  ratios  of  1.6  (166 
Hz/99  Hz)  and  1.7  (166  Hz/93  Hz),  making  it  unlikely  that  interaction  of  the  fundamental 
frequencies  created  strong  distortion  products  in  the  cochlea  (Dhar  et  al.,  2009,  2005;  Knight  and 
Kemp,  2001).  Moreover,  the  responses  of  trained  musicians  were  significantly  enhanced 
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compared  with  those  of  novice  listeners,  implying  experience-  based  differences  that  would  not 
have  arisen  at  the  level  of  the  cochlea  or  auditory  nerve  (e.g.,  Large,  Kozloski,  &  Crawford,  1998; 
Laudanski  et  al.,  2010). 

The  stimuli  from  the  Lee  et  al.  (2009)  study  were  used  as  input  to  a  cochlear  model  (see 
Section  5.1),  which  in  turn  provided  input  to  two  brainstem  network  layers.  The  characteristic 
frequencies  of  the  cochlear  layer  and  both  brainstem  layers  spanned  four  octaves  with  99 
oscillators  per  octave.  Thus,  each  layer  included  397  oscillators,  with  characteristic  frequencies 
ranging  from  64  Hz  to  1024  Hz,  encompassing  the  range  of  frequencies  for  which  time-locked 
responses  have  been  observed  in  midbrain  physiology  (Langner,  1992).  The  cochlear  model 
includes  a  middle  ear  filter  and  simulates  the  basilar  membrane  and  the  organ  of  Corti  (cf. 
Jlilicher  et  al.,  2001  B).  The  cochlea  is  connected  to  the  first  brainstem  layer,  representing  the 
cochlear  nucleus  (CN),  and  the  CN  is  connected  to  the  second  brainstem  layer,  representing  the 
inferior  colliculus/  lateral  lemniscus  (IC/LL). 

We  modeled  FFRs  from  Lee  et  al.,  responses  to  both  the  consonant  and  dissonant  intervals. 
The  stimulus  was  input  to  the  cochlea,  all  oscillator  equations  were  numerically  integrated  for 
the  length  of  the  stimulus,  and  the  responses  in  all  layers  were  stored.  To  compute  the  model 
brainstem  FFR,  the  responses  of  all  oscillators  in  each  layer  were  averaged,  leaving  a  single  time 
series  for  each  layer.  The  model  FFR  was  a  weighted  average  of  the  cochlea,  CN,  and  IC/LL 


Figure  5.2.  Comparisons  of  model  predictions  and  auditory  brainstem  responses  of  nonmusicians  to  (A)  the 
consonant  interval  (99  Hz,  166  Hz)  and  (B)  the  dissonant  interval  (93  Hz,  166  HzThe  labels  above  each  spectral 
component  refer  only  to  their  specific  frequencies  as  functions  of  the  primaries,  and  do  not  necessarily  reflect  the 
aeneratina  processes  of  those  components. 

layers.  This  weighted  average  was  filtered  (3rd-order  digital  Butterworth  low-pass,  450  Hz  cutoff) 
to  account  for  the  lowpass  effect  of  the  skull,  meninges,  and  scalp  on  the  FFR.  Finally,  the 
resultant  model  time  series  were  averaged  and  fast  Fourier  transformed  to  produce  a  model  fit. 
The  model  fits,  for  the  consonant  and  dissonant  intervals  were  optimized  with  a  single  degree  of 
freedom  through  a  series  of  simulations.  In  the  simulations,  the  parameter  s  (Equation  1)  was 
systematically  varied  between  zero  and  one  to  yield  the  highest  correlation  for  each  fit. 

Fig.  5.2  shows  the  predicted  brainstem  responses.  The  canonical  population  response 
predicts  each  peak  in  the  brainstem  response  with  remarkable  accuracy  for  both  intervals 
(consonant:  R2=0.77,  p<0.0001 ;  dissonant:  R2=0.67,  p<0.0001 ).  The  more  subtle  differences 
could  not  be  accounted  for  by  manipulation  of  s  alone,  implying  that  incorporation  of  other 
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network  properties  into  the  model,  i.e.,  synaptic  coupling,  will  be  necessary  to  explain  the 
responses  of  listeners.  This  is  consistent  with  the  interpretation  that  the  refinement  of  auditory 
sensory  encoding  is  driven  by  synaptic  plasticity  that  links  learned  representations  to  the  neural 
encoding  of  acoustic  features  (Lee  et  ah,  2009). 

Summary:  A  single  parameter  of  a  model  brainstem  network  was  manipulated  to  fit  Lee  et 
al.’s  (2009)  brainstem  FFR  data.  The  parsimony  of  the  model,  its  basis  in  neurophysiological 
observations  of  mode-locking,  and  the  quality  of  the  fits  all  speak  to  the  potential  of  this 
theoretical  approach.  We  are  currently  developing  improvements  to  enable  more  comprehensive 
simulations  of  the  early  auditory  system  that  include  all  relevant  aspects  of  cochlear  dynamics,  as 
well  as  parameterization  of  CN  and  IC/LL  dynamics.  In  comprehensive  models,  parameter  fitting 
is  a  significant  issue.  We  are  also  exploring  models  of  synaptic  plasticity  for  neural  oscillator 
networks  (see  Equation  2  and  Section  4)  in  order  to  better  explain  the  responses  of  trained 
listeners.  This  approach  may  lead  to  an  understanding  of  general  neural  signal  processing 
principles  underlying  music  and  pitch  perception.  Moreover,  canonical  analysis  of  plasticity  in 
neural  oscillator  networks  may  help  us  to  understand  the  role  of  learning  in  modulating  these 
responses. 

5.3  Cortical  Synchronization  to  Complex  Rhythms 

We  have  developed  a  model  that  makes  detailed,  quantitative  predictions  about  cortical 
population  rhythms  hypothesized  to  underlie  human  rhythm  perception.  We  directly  tested  the 
predictions  in  behavioral  and  neurophysiological  experiments.  In  particular,  our  model  makes 
predictions  about  the  perception  of  pulse  and  meter  in  complex  rhythms  as  well  as  activity  within 
auditory  and  motor  systems. 

Pulse  and  meter  perception  arise  from  complex  interactions  within  a  widespread 
auditory-motor  network  (Lee,  Skoe,  Kraus,  &  Ashley,  2009).  EEG  and  MEG  studies  have 
directly  tested  the  neural  dynamics  that  emerge  during  rhythm  perception.  One  group  of  studies 
has  looked  at  how  steady  state  evoked  potentials  (SS-EP)  are  affected  by  pulse  perceptions  and 
by  different  tempi.  Will  and  Berg  (e.g.,  Chen,  Penhune,  &  Zatorre,  2008;  Grahn  &  Rowe,  2009) 
reported  substantial  SS-EP  responses  to  isochronous  stimuli,  but  only  in  the  delta/theta  range, 
with  the  strongest  SS-EP  response  for  2Hz  stimuli,  a  rate  which  corresponds  to  the  optimal 
pulse-tempo  identified  in  numerous  listening,  tapping  synchronization,  and  event-interval 
discrimination  experiments  (2007).  Another  study  found  that  while  periodic  rhythm  elicited  a 
sustained  response  at  the  rate  of  the  stimulus,  meter  imagery  elicited  an  additional  subharmonic 
resonance  corresponding  to  the  metric  structural  interpretation  (McAuley  &  Jones,  2003). 
Similarly,  the  amplitude  of  the  SS-EPs  at  pulse  and  meter  frequencies  of  complex  rhythms  is 
selectively  enhanced  (S.  Nozaradan  et  al.,  2011).  While  these  studies  support  the  notion  that 
neural  activity  synchronizes  to  rhythms,  synchronization  cannot  be  unambiguously  attributed  to 
oscillatory  entrainment  because  they  do  not  rule  out  the  possibility  that  SS-EPs  may  be  driven 
bottom-up  by  the  stimulus  itself. 

Other  EEG/MEG  studies  have  demonstrated  synchronization  of  high-frequency  power 
modulations  (in  the  beta/gamma-band  range  15-50Hz)  to  the  temporal  structure  of  sounds.  One 
EEG  study  found  that  fluctuations  in  induced  beta-  and  gamma-band  power  synchronized  with 
periodic  and  metrical  rhythms,  and  was  unaltered  even  when  sounds  were  omitted,  emphasizing 
its  top-down  anticipatory  (rather  than  bottom-up  reflexive)  nature  (Sylvie  Nozaradan,  2012). 
Similarly,  MEG  studies  found  anticipatory  induced  beta-band  responses  for  periodic  and  metrical 
sequences,  but  not  for  randomly  timed  sequences  (Snyder  &  Large,  2005),  and  subharmonic 
responses  were  observed  in  induced  beta  band  activity  when  subjects  were  instructed  to  impose  a 
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subjective  meter  on  a  periodic  stimulus,  closely  resembling  those  produced  by  physical  accents 
(Fujioka,  Large,  Trainor,  &  Ross,  2009).  MEG  source  analysis  has  revealed  beta-band 
interactions  in  auditory  and  motor  networks  during  musical  rhythm  processing  (Iversen,  Repp,  & 
Patel,  2009).  Given  the  role  of  beta  in  motor  processing  and  long-range  intra-cortical  interaction, 
these  findings  are  consistent  with  the  idea  that  the  motor  system  influences  the  perception  of 
sound,  even  in  the  absence  of  overt  movement.  Thus,  beta  and  gamma  band  responses  to 
auditory  rhythms  as  well  as  SS-EPs  are  in  line  with  the  basic  prediction  that  neural  rhythms 
synchronize  to  the  pulse. 

While  empirical  observations  demonstrate  neural  synchrony,  the  hypothesis  of  entrained 
neuronal  oscillation  remains  controversial.  The  most  straightforward  objection  is  that  1)  apparent 
entrainment  of  SS-EPs  reflects  overlapping  of  transient  responses  such  as  the  N1-P2  complex 
(Sussman,  Steinschneider,  Gumenyuk,  Grushko,  &  Lawson,  2008;  Tremblay,  Billings,  &  Rohila, 
2004).  A  related  possibility  is  that  2)  observed  synchrony  may  be  a  passive  response,  like  that  of 
a  linear  band-pass  filter.  If  synchrony  is  passive  then  3)  observed  neural  processes  would  not  be 
capable  of  the  cognitive  computations  necessary  for  structure  perception  (Patel  &  Iversen,  2014). 
Fortunately,  it  is  possible  to  identify  active  responses  because  entrainment  reflects  interaction  of 
the  stimulus  with  intrinsic  neural  dynamics. 

5.3.1  The  Auditory-Motor  Resonance  Model 


We  used  two  gradient  frequency  networks  (see  Equation  1)  to  model  the  functional 
coupling  of  auditory-motor  networks  observed  in  rhythm  perception  tasks  without  a  motor 
component.  The  sensory  network  is  intended  to  capture  auditory  cortical  entrainment,  while  the 
motor  network  is  intended  to  capture  the 
dynamics  of  a  broadly  distributed 
network  including  basal  ganglia  and 
cortical  areas.  The  sensory  network  takes 
a  rhythmic  input,  sends  output  to  a  motor 
network,  and  the  motor  network  send 
input  back  to  the  sensory  network 
(Large  et  al.,  2015).  Connections  within 
and  between  networks  are  assumed  to  be 
plastic  and  tuned  by  musical 
enculturation  (cf.,  Hannon  and  Trehub, 

2005;  see  also  Section  4). 

Although  the  model  makes  only 
general  assumptions  regarding 
underlying  neural  structures  (e.g.,  Chen 
et  al.,  2008),  it  makes  strong 
commitments  about  the  oscillatory 
dynamics  of  auditory-motor  interactions 
(Will  and  Berg,  2007;  Fujioka  et  ah, 

2012;  Nozaradan  et  ah,  2013).  The 
sensory  oscillators  are  tuned  to  operate 
near  a  Hopf  bifurcation;  the  motor 
oscillators  are  tuned  to  operate  near  a 
double  limit  cycle  bifurcation  (see 
Section  2).  The  double  limit  cycle 


frequency 


Figure  5.3.  The  network  was  stimulated  with  (a)  an  isochronous 
rhythm  and  (b)  a  “missing  pulse”  rhythm.  Output  of  the  sensory 
network  is  in  green  and  output  of  the  sensory  network  is  in  blue. 
The  mean  field  time  series  (left)  was  obtained  by  summing  the 
output  of  all  the  oscillators  in  the  network  over  time.  The  SS-EP 
was  obtained  by  Fourier  analysis  (FFT)  of  the  mean  field  time 
series.  The  Fourier  analysis  (FFT)  of  the  stimulus  envelope  is 
shown  in  black  on  the  sensory  SS-EP  axis. 
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regime  of  the  motor  network  means  that  the  model  can  capture  synchronization-continuation 
behavior,  continuing  to  produce  rhythmic  behavior  after  the  stimulus.  To  predict  mean  field  time 
series  as  observed  in  EEG  recordings  (e.g..  Will  and  Berg,  2007;  Stefanics  et  al.,  2010),  we  sum 
the  output  of  all  oscillators  in  each  network  (Figure  5.3,  left).  To  predict  steady  state  evoked 
potentials  (SS-EPs)  we  take  a  frequency  analysis  (DFT)  of  the  mean  field  (Figure  5.3,  right). 

As  shown  in  Figure  5. 3 A,  for  a  periodic  stimulus,  both  sensory  and  motor  networks 
produce  synchronized  oscillations  at  the  pulse  frequency,  and  generate  harmonics  (Repp,  2008) 
and  subharmonics  (Vos,  1973;  Bolton,  1894;  Nozaradan  et  al.,  2011).  In  the  case  of  a  complex 
rhythm,  however,  it  becomes  clear  that  the  two  networks  are  doing  something  quite  different 
from  one  another.  The  mean  field  time  trace  for  the  sensory  network  represents  the  input  rhythm 
rather  faithfully,  producing  well-defined  pulses  at  input  event  times.  By  contrast,  the  motor 
network  entrains  at  the  pulse  frequency.  The  rhythm  itself  contains  no  energy  at  the  pulse 
frequency  (or  its  second  subharmonic;  DFT  in  Figure  5.3B,  SS-EP,  solid  black),  however,  in  the 
motor  network  the  strongest  response  is  found  at  the  pulse  frequency.  In  other  words,  the 
development  of  the  pulse  percept 
depends  on  the  interaction  of  these 
two  oscillatory  systems. 

This  predicts  that  an 
oscillatory  network  interaction  can 
lead  to  spontaneous  pulse  induction 
in  complex  rhythms — even  in  the 
most  extreme  case  of  a  rhythm  for 
which  there  is  no  energy  at  the  pulse 
frequency.  Thus,  the  theoretical 
prediction  is  that  pulse  may  be 
perceived  at  a  frequency  that  is  not 
physically  present  in  the  rhythmic 
stimulus  (Large,  2010;  Velasco  and 
Large,  2011). 

5.3.2  Behavioral  Data 

To  test  behavioral 
predictions,  we  asked  participants  to 
listen  to  eleven  rhythms  ranging  from 
isochronous  to  highly  complex. 

Participants  were  instructed  to  listen 
to  each  rhythm  until  they  heard  a 
steady  pulse,  and  then  tap  along  with 
the  rhythm  at  that  rate.  Rhythms  were 
presented  at  five  levels  of  complexity 
(0-4)  and  at  five  different  tempi  (i.e., 
pulse  frequencies):  2.28Hz  (420  ms), 

2.17  Hz  (460  ms),  2  Hz  (500  ms), 

1.85  Hz  (540  ms),  and  1.72  Hz  (580 
ms).  The  simplest  rhythms 
(complexity  0)  were  isochronous,  the 
most  complex  rhythms  contained  no 
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Figure  5.4.  (a)  Tapping  frequencies  were  normalized  to  2  Hz  to 
allow  comparison  between  trials  at  different  tempos.  Tapping 
frequency  distributions  (red  histograms)  were  computed  by 
binning  normalized  instantaneous  tapping  frequencies  from  0 
Hz  to  5.00  Hz  in  bins  widths  of  0.05  Hz.  Distributions  were 
computed  for  each  rhythm  separately,  including  every  tap 
interval  across  trials.  Black  lines  show  amplitude  spectrum  of 
the  stimulus  envelope  for  comparison,  (b)  Circular  means  of  tap 
phases  for  each  trial  (blue  circles)  and  grand  mean  for  each 
comolexitv  level  (red  line). 
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spectral  energy  at  the  pulse  frequency  (complexity  4).  Different  combinations  of  tempo  and 
rhythm  were  presented  in  a  pseudorandom  order  such  that  consecutive  trials  always  had  different 
rhythms,  and  different  tempos.  Thus,  participants  were  forced  to  find  both  the  frequency  and  the 
phase  of  the  pulse  anew  for  each  rhythm;  they  could  not  simply  tap  at  the  same  tempo 
throughout  the  experiment. 

We  measured  instantaneous  tapping  frequency  distributions  to  determine  whether 
subjects  induced  a  pulse  at  the  intended  frequency.  Instantaneous  tapping  frequency  was 
computed  as  1/ITI  (ITI  =  inter-tap  interval  in  seconds)  and  tapping  frequencies  were  normalized 
to  a  frequency  of  2  Hz  so  they  could  be  combined  into  a  single  distribution  at  each  level  of 
complexity.  Spectral  analysis  (DFT)  of  the  stimulus  rhythms  (Figure  5.4A,  black)  shows  that  at 
the  hypothetical  pulse  frequency  amplitude  decreases  with  increasing  complexity.  At  complexity 
level  4,  the  amplitude  is  precisely  zero  at  2  and  1  Hz  for  each  rhythm.  Normalized  instantaneous 
tapping  frequency  (Figure  5. 4 A,  red  histogram)  displays  a  main  peak  at  the  normalized  pulse 
frequency  of  2Hz  for  all  rhythms  at  all  levels  of  complexity,  with  lesser  peaks  at  1  and  0.5Hz, 
and  for  some  rhythms,  a  diffuse  peak  around  4Hz.  Thus,  the  participants  most  often  tapped  the 
predicted  pulse  frequency  even  for  the  most  complex  rhythms,  which  had  no  spectral  amplitude 
at  that  frequency. 

Next,  we  examined  synchronization  for  each  trial.  The  sequence  of  tap  times  was 
converted  into  a  sequence  of  phases  relative  to  the  predicted  pulse  frequency,  and  the  circular 
mean  was  computed  for  each  trial  (Batschelet,  1981;  Figure  5.4B;  blue  circles).  The  grand  mean 
was  then  computed  for  each  complexity  level  (Figure  5.4B,  red  line).  As  predicted  by  the  model, 
participants  synchronized — either  in-phase  or  anti-phase — predominantly  at  the  missing  pulse 
frequency.  This  behavior  is  consistent  with  the  prediction  that  formation  of  the  pulse  percept 
arises  due  to  entrainment  of  emergent  neuronal  oscillations.  It  also  rules  out  the  potential 
alternatives  that  synchronization  is  merely  a  consequence  of  a  common  rhythmic  input,  or  that 
the  pulse  percept  may  arise  due  to  linear  resonance.  Theoretically  speaking,  it  is  critical  to 
distinguish  the  role  of  a  common  stimulus  frequency  from  the  intrinsic  dynamics  of  an  emergent 
oscillation  (Whittington  et  al.,  2000),  and  the  missing  pulse  rhythms  used  here  enabled  us  to 
dissociate  the  two. 


5.3.3  Neural  Data 

The  theoretical  model  also  makes  strong 
physiological  predictions  that  must  also  be  put  to 
empirical  test.  The  prediction  is  that  in  listening  to 
these  ‘missing  pulse’  rhythms,  the  pulse  frequency 
will  be  observed  in  auditory  and  motor  activity 
despite  its  absence  from  the  stimulus  acoustics.  Fig 
5.5  shows  MEG  data  from  n=8  participants, 
listening  passively  to  two  of  the  complex  ‘missing 
pulse’  stimuli  (MP1  and  MP2)  used  in  the  model 
simulations  (Tal  et  al,  submitted).  The  spectra  of 
responses  from  right  and  left  auditory  cortices  show 
the  emergence  of  peaks  at  the  2Hz  ‘missing  pulse’ 
(and  in  one  case  its  1Hz  subharmonic;  p<0.05  for 
all  relative  to  adjacent  frequencies;  regions  of 
interest  selected  based  on  an  independent  auditory 


Frequency  (Hz) 

Figure  5.5.  (a)  Synchronized  neural  activity  was 
found  in  left  and  right  auditory  cortices  in 
response  to  both  complex  rhythms,  (b)  Energy 
was  found  at  the  missing  pulse  (here,  2  Hz)  and 
its  second  subharmonic  (1  Hz)  consistent 
neither  of  which  were  present  in  the  stimuli. 
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localizer  task;  Source  localization  performed  using  SAM  beamforming).  Moreover,  we  found 
that  the  magnitude  of  the  ‘missing  pulse’  response  is  correlated  with  pulse  perception  across 
individuals,  such  that  participants  who  perceived  the  pulse  quickly  had  stronger  responses  than 
participants  who  took  longer  to  perceive  a  pulse  (MP1:  r=-0.81,  p<0.005;  MP2:  r=-0.6,  p<0.05; 
data  not  shown).  These  findings  confirm  the  falsifiable  physiological  predictions  of  the  model, 
they  show  that  neural  frequencies  correspond  to  listener’s  perceptions,  and  they  conclusively  rule 
out  transient  or  passive  explanations  of  neural  synchrony.  Thus,  this  experiment  provided  strong 
evidence  for  the  viability  of  the  model  as  a  theoretical  framework  for  explaining  rhythm 
perception. 

Summary:  We  have  made  significant  progress  in  understanding  the  role  of  neural 
oscillations  and  the  neural  structures  that  support  synchronized  responses  to  musical  rhythm.  Our 
neurodynamic  model  that  shows  how  self-organization  of  oscillations  in  interacting  sensory  and 
motor  networks  could  be  responsible  for  the  formation  of  the  pulse  percept  in  complex  rhythms. 
In  a  pulse  synchronization  study,  we  tested  the  model’s  key  prediction  that  pulse  can  be 
perceived  at  a  frequency  for  which  no  spectral  energy  is  present  in  the  amplitude  envelope  of  the 
acoustic  rhythm.  The  result  shows  that  participants  perceive  the  pulse  at  the  theoretically 
predicted  frequency.  Synchronized  neural  activity  consistent  with  model  predictions  has  also 
been  observed  in  auditory  cortex,  providing  strong  evidence  for  the  viability  of  this  theoretical 
framework  for  explaining  rhythm  perception.  Our  model  is  one  of  the  few  consistent  with 
neurophysiological  evidence  on  the  role  of  neural  oscillation,  and  it  explains  a  phenomenon  that 
other  computational  models  fail  to  explain.  Because  it  is  based  on  a  canonical  model,  the 
predictions  hold  for  an  entire  family  of  dynamical  systems,  not  only  a  specific  one.  Thus,  this 
model  provides  a  theoretical  link  between  oscillatory  neurodynamics  and  the  induction  of  pulse 
and  meter  in  musical  rhythm. 

6.  Concluding  Remarks 

In  this  project,  we  developed  a  theoretical  framework  for  auditory  neural  processing  and 
auditory  perception.  We  modeled  the  auditory  system  as  a  dynamical  system  consisting  of 
oscillatory  networks,  and  auditory  perception  as  stable  dynamic  patterns  formed  in  the  networks 
in  response  to  acoustic  signals.  We  developed  GrFNNs,  generic  models  that  capture  the 
neurocomputational  properties  of  a  family  of  neurophysiological  models  using  bifurcation  theory. 
We  made  significant  progress  in  understanding  the  signal  processing,  pattern  formation  and 
plasticity  in  GrFNNs.  We  developed  three  models  that  exploit  these  properties  to  model 
important  aspects  of  auditory  neurophysiology  and  auditory  perception.  Future  modeling  efforts 
based  on  canonical  dynamical  systems  could  bring  us  closer  to  understanding  fundamental 
mechanisms  of  hearing,  communication,  and  auditory  system  development. 

In  addition  to  these  accomplishments,  we  produced  a  computational  framework  for  GrFNNs 
in  Matlab,  sped  up  our  computational  simulations  using  GPU  acceleration,  and  created  a  C++ 
version  of  the  GrFNN  code  to  develop  end-user  applications  that  run  on  CPUs,  GPU,  mobile 
platforms  and  embedded  devices.  We  anticipate  that  this  computational  framework  will  power 
the  next  generation  of  auditory  processing  hardware  and  software.  In  ongoing  work,  we  are 
exploring  the  abilities  of  such  networks  to  address  the  auditory  scene  analysis  problem. 
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